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Abstract 

Determining  the  fate  and  transport  of  JP-8  jet  fuel  is  a  complex  and  important 
problem.  As  part  of  the  startup  procedures  for  jet  engines,  fuel  is  passed  through 
aircraft  engines  before  combustion  is  initiated.  Because  of  the  extremely  low  tem¬ 
peratures  at  northern  tier  Air  Force  bases,  the  unburned  fuel  does  not  evaporate 
readily  and  may  come  into  contact  with  ground  crew.  To  determine  the  amount 
and  duration  of  contaminant  contact,  the  evaporation  of  the  emitted  fuel  must  be 
modeled.  The  amount  and  composition  of  the  fuel  upon  reaching  the  ground  crew 
may  be  determined  by  droplet  evaporation  models  that  have  already  been  developed. 
The  evaporation  of  the  fuel  after  adhering  to  the  skin  needs  to  be  modeled.  This 
knowledge  of  the  fuel’s  fate  may  then  be  used  to  determine  source  terms  for  use  in 
toxicological  studies. 

This  research  involves  the  comparison  of  two  existing  droplet  evaporation  mod¬ 
els  and  the  calculation  of  the  evaporation  of  a  film  of  jet  fuel  from  a  surface.  The 
existing  models  are  compared  in  order  to  make  recommendations  on  which  model  to 
use  to  predict  the  amount  and  composition  of  fuel  reaching  the  ground  crew.  To  make 
the  surface  evaporation  problem  amenable  to  modeling,  simplifying  assumptions  are 
made.  The  fuel  is  assumed  to  be  a  uniformly  distributed  mixture  of  representative 
hydrocarbon  groups.  Due  to  the  complexity  of  the  mixture  of  aviation  fuels,  a  mix¬ 
ture  of  the  predominant  species  were  chosen  as  representatives  to  approximate  the 
physical  behavior  of  the  entire  fuel  mixture. 

The  goal  of  this  research  is  to  determine  the  most  appropriate  model  for  pre¬ 
dicting  the  amount  and  composition  of  jet  fuel  reaching  the  ground  crew  and  to 
extend  the  more  appropriate  fuel  droplet  evaporation  model  to  describe  the  evapo¬ 
ration  of  a  film  of  fuel  from  a  surface.  A  validation  of  the  resultant  model  is  then 
performed  by  comparing  the  calculations  to  experimental  data. 
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Evaporation  of  Jet  Fuels 


/.  Introduction 

1.1  Overview 

Determining  the  fate  and  transport  of  JP-8  jet  fuel  is  a  complex  and  important 
problem.  Normal  air  force  operations  sometimes  require  fuel  to  be  released  into  the 
atmosphere.  Fuel  is  purposely  jettisoned  from  aircraft  during  emergency  situations 
to  reduce  the  risk  of  fire  or  to  make  it  easier  to  land  safely.  Fuel  is  also  passed 
through  aircraft  engines  without  combustion  in  order  to  prepare  for  flight. 

Previous  research  into  the  evaporation  and  settling  to  the  ground  of  jettisoned 
fuel  has  led  to  the  development  of  two  droplet  evaporation  and  transport  models  to 
determine  the  amount  and  concentration  of  contaminant  that  reaches  the  ground 
after  jettisoning.  This  knowledge  of  the  fuel’s  fate  has  been  used  to  determine 
minimum  altitudes  for  jettisoning  without  contamination  of  the  ground. 

More  recently,  there  have  been  reports  of  ground  crews  being  exposed  to  un¬ 
burned  fuel  during  the  cold  start  of  jet  engines  at  northern  tier  Air  Force  bases. 
Fuel  is  passed  through  the  aircraft  engines  to  prepare  for  flight.  When  the  ambient 
temperature  is  extremely  low,  as  it  is  at  northern  tier  Air  Force  bases,  the  fuel  does 
not  evaporate  quickly.  The  amount  and  composition  of  fuel  reaching  the  ground 
crew  can  be  predicted  using  existing  droplet  evaporation  models  (3),  (16).  However, 
there  is  a  requirement  to  model  the  evaporation  of  the  fuel  from  the  ground  crew’s 
skin  and  clothing.  The  evaporation  data  could  then  be  used  to  determine  source 
terms  for  use  in  toxicological  and  pharmacokinetic  studies  such  as  determining  the 
amount  of  fuel  that  is  absorbed  through  the  skin.  Integrating  models  of  evaporation 
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and  absorption  of  the  fuel  through  skin  may  then  provide  a  quantification  of  the 
extent  of  exposure  to  the  ground  crew. 

1.2  Problem 

In  order  to  assess  the  exposure  of  the  fuel  on  the  contacted  ground  crew,  evap¬ 
oration  calculations  are  needed.  It  is  unclear  whether  or  not  the  several  droplet 
evaporation  algorithms  that  have  been  developed  produce  equivalent  results.  Com¬ 
parisons  must  be  performed  to  determine  if  there  is  a  significant  difference  in  the 
calculations  and  to  determine  which  model  should  be  used  to  predict  the  amount 
and  composition  of  fuel  reaching  the  ground  crew.  The  ease  of  application  of  these 
models  to  newly  developed  fuels  also  needs  to  be  considered  as  the  fuel  used  by  the 
Air  Force  is  continually  changing. 

The  fuel  will  not  likely  consist  of  spherical  droplets  after  adhering  to  an  ob¬ 
ject.  Therefore,  droplet  evaporation  studies  that  assume  spherical  droplets  are  not 
applicable  to  the  fuel  after  it  reaches  the  crew. 

1.3  Scope 

The  scope  of  this  research  involves  the  numerical  computation  of  the  evapora¬ 
tion  of  a  film  of  jet  fuel  from  a  surface  and  the  comparison  of  existing  fuel  droplet 
evaporation  descriptions.  The  droplet  evaporation  calculations  will  be  examined  in 
an  attempt  to  determine  which  model  should  be  used  to  predict  the  amount  and 
composition  of  the  fuel  that  reaches  the  ground  crew.  The  surface  evaporation  cal¬ 
culation  will  be  based  on  the  previous  studies  of  droplet  evaporation.  The  droplet 
evaporation  models  will  be  modified  as  necessary  to  describe  the  evaporation  of  the 
fuel  from  a  flat  surface. 

The  fuel  is  assumed  to  be  a  uniformly  distributed  mixture  of  representative 
hydrocarbon  groups.  Due  to  the  complexity  of  the  mixture  of  aviation  fuels,  a 
mixture  of  the  predominant  species  were  chosen  as  representatives  to  approximate 
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the  physical  behavior  of  the  entire  fuel  mixture.  The  evaporation  of  the  fuel  is 
complicated  by  the  fact  that  the  composition  of  the  droplet  is  continually  changing 
because  the  more  volatile  components  evaporate  more  quickly  than  the  less  volatile 
components.  Furthermore,  the  temperature  of  the  fuel  is  continually  changing  which 
affects  the  evaporation  rate.  The  fuel  temperature  approaches  equilibrium  with  the 
ambient  air  temperature  and  is  simultaneously  cooled  through  evaporative  cooling. 
These  complications  lead  to  a  nonlinear  system  of  equations  describing  evaporation 
that  must  be  solved  numerically. 

1.4  Summary  of  Thesis 

The  thesis  is  organized  as  follows:  Chapter  1,  Introduction;  Chapter  2,  Background 
of  evaporation  theory;  Chapter  3,  Comparison  of  droplet  evaporation  algorithms; 
Chapter  4,  Surface  evaporation  algorithm;  Chapter  5,  Conclusions  and  Recommen¬ 
dations;  Bibliography. 
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IL  Background 


2.1  Overview 

Evaporation  is  a  heat  and  mass  transport  phenomenon.  Thus,  in  order  to  un¬ 
derstand  and  predict  rates  of  evaporation,  one  must  understand  the  rates  associated 
with  these  processes.  In  the  following  sections,  heat  and  mass  transfer  equations  that 
describe  evaporation  of  a  liquid  mixture  are  developed.  The  equations  presented  in 
the  droplet  evaporation  descriptions  begin  with  expressions  for  the  time  derivatives 
of  droplet  or  species  mass,  temperature,  and  radius.  The  heat  and  mass  transfer 
equations  are  presented  here  in  order  to  illustrate  that  the  calculations  of  droplet 
evaporation  are  derived  from  first  principles.  Following  the  development  in  section 
2.2,  the  droplet  evaporation  calculations  are  presented  and  are  referenced  to  the 
equations  in  section  2.2. 

The  principal  references  used  for  the  heat  and  mass  transfer  rates  presented 
in  section  2.2  are  Bird,  Stewart,  and  Lightfoot(l)  and  Incropera  and  De  Witt(7). 
Supplemental  material  is  also  taken  from  Clark  (2)  and  Fuchs  (5).  The  principal 
references  for  the  theory  applied  in  the  combustion  model  are  Classman  (6),  Kuo(8), 
and  Williams(18). 

2.2  Physics  of  Evaporation 

2.2.1  Mass  Transfer. 

2.2. 1.1  Definition  of  Terms.  We  begin  by  defining  the  concentra¬ 
tions,  velocities,  and  fluxes  that  will  be  used  to  develop  the  mass  transfer  equation. 

The  fuel  is  composed  of  several  components  that  we  refer  to  as  species.  The 
concentration  of  species  i  may  be  expressed  in  several  ways.  The  mass  concentration, 
pi,  is  the  mass  m  of  species  i  per  unit  volume  of  solution  V.  The  molar  concentration, 
Ci  =  Pi/ Wi,  where  Wi  is  the  molecular  weight  of  species  i,  is  the  number  of  moles  per 
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unit  volume  of  solution.  The  mass  fraction,  y*  =  pi/ p,  is  the  mass  concentration  of  i 
divided  by  the  density  of  the  solution.  The  mole  fraction,  Xi  =  Cj/c,  is  calculated  as 
the  molar  concentration  of  i  divided  by  the  total  number  of  moles  per  unit  volume 
of  solution,  c. 

Likewise,  the  velocity  of  a  species  moving  through  a  solution  may  be  expressed 
in  different  ways.  Here,  it  is  defined  in  terms  of  the  mass  because  we  want  an 
expression  for  the  mass  flux.  The  velocity  refers  to  the  sum  of  the  velocities  of  all 
molecules  of  a  species  divided  by  the  total  number  of  the  molecules  in  a  small  volume. 
The  mass  average  velocity,  v,  of  a  mixture  of  N  species  is  defined  by; 

V  =  (2.1) 

P 

where  p  =  Pi  is  the  total  mass  concentration. 

The  mass  flux  of  species  A  relative  to  stationary  coordinates  is  defined  as: 

Ha  =  PaVa,  (2.2) 

The  mass  flux  of  A  relative  to  coordinates  moving  at  the  average  mass  velocity  of 


the  mixture  is  defined  as: 


=  Pa{^a  -  v) 


Replacing  Pa^a  in  the  definition  of  Ja  with  ua  results  in: 


jA  =  nA  -  Pav 


Substituting  the  definition  of  v  yields: 


Ja  =  “A  —  Pa 


Eili  Pi^i 
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For  a  binary  mixture  composed  of  species  A  and  B,  we  may  use  the  definition  of  the 
mass  fraction  of  species  li  =  Pi/p,  and  the  definition  the  mass  fiux  of  species  i, 
=  piVi,  to  replace  equation  (2.5)  with: 


-  YA{nA  +  ns)  (2.6) 

2.2. 1.2  Mass  Diffusion.  To  describe  the  diffusion  of  species  A  through 
the  mixture,  Pick’s  first  law  is  applied  to  give  the  mass  flux  of  A  with  respect  to  the 
mass  average  velocity: 

=  -pVVYa  (2.7) 

where  V  is  the  diffusion  coefficient  and  VYa  is  the  concentration  gradient  of  A. 
Notice  that  the  sign  of  the  right  hand  side  of  equation  (2.7)  is  negative.  That  means 
that  the  mass  will  fiow  in  the  direction  of  steepest  concentration  descent.  That  is, 
A  will  move  from  an  area  of  high  concentration  to  an  area  of  low  concentration. 

Rearranging  equation  (2.6)  and  replacing  Ja  using  Pick’s  law  yields: 


Ha  =  YAivkA  +  ub)  -  pWYa  (2.8) 

where  pWYa  is  the  contribution  to  the  mass  flux  by  ordinary  diffusion  and  FA(nA  + 
iib)  is  the  contribution  due  to  the  diflPusion  induced  bulk  movement  of  A  and  B.  To 
see  that  the  bulk  contribution  is  physically  sensible,  consider  what  happens  as  A 
begins  to  evaporate  and  diffuse  away  from  the  surface.  There  is  clearly  a  flux  of  A 
away  from  the  liquid.  If  we  assume  that  gas  B  is  insoluble  in  liquid  A,  there  also 
must  be  a  flux  of  B  away  from  the  liquid  as  it  moves  to  make  room  for  the  additional 
molecules  of  A.  The  result  is  a  net  flux  of  both  A  and  B  away  from  the  surface. 

Pick’s  law  describes  diffusion  when  there  are  no  convective  affects.  This  has 
application  to  problems  where  the  medium  through  which  A  is  diffusing  is  motionless. 
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To  handle  problems  where  the  medium  is  in  motion,  we  must  consider  mass  transport 
by  diffusion  and  convection. 


2.2. 1.3  Mass  Diffusion-Convection.  Consider  mass  transfer  across  a 
liquid- vapor  phase  boundary  as  in  Figure  2.1  where  liquid  A  is  evaporating  into  a 
gas  B,  and  gas  B  is  moving  with  a  constant  velocity  u. 


Figure  2.1  Depiction  of  evaporation  process. 


The  motion  of  B  across  the  surface  causes  convective  mass  transfer  to  occur 
in  addition  to  the  diffusion  mass  transfer.  Problems  of  this  sort  are  generally  too 
complex  to  solve  analytically.  Instead,  a  mass  transfer  coefficient,  k,  is  defined  as  in 
Bird(l)  in  terms  of  the  rate  of  diffusion  normal  to  the  gas-liquid  interface  as  follows: 

jAo  =  kAYA  (2.9) 

where  jao  is  the  flux  of  gas  A  into  the  gas  B  at  the  liquid- vapor  boundary,  with  the 
zero  subscript  to  denote  conditions  at  the  liquid-vapor  boundary,  and  Al^i  is  the 
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difference  between  the  gas  phase  concentration  of  A  at  the  liquid-vapor  boundary 
and  in  the  approaching  stream  of  gas  B.  Here,  k  has  been  averaged  over  the  surface 
area  of  the  liquid  A.  Notice  that  the  sign  of  the  right  hand  side  of  equation  (2.9)  is 
positive  because  jao  is  the  mass  evaporating  per  unit  area  per  unit  time. 

Evaluating  equation  (2.6)  at  the  surface  of  the  liquid,  denoting  this  with  a  zero 
subscript,  and  rearranging  results  in: 

^Ao  =  yAai^AQ  +  nj5o)  +  3  Ad 
Substituting  jAd  from  equation  (2.9)  yields: 

“^Ad  =  yAdip>Ad  +  f^Bd)  +  kAYA  (2.10) 

where  n^o  is  the  mass  flux  at  the  liquid  surface  with  respect  to  fixed  coordinates. 

The  mass  transfer  coefficient  is  an  empirically  derived  parameter  that  repre¬ 
sents  all  convective  effects.  At  high  mass  transfer  rates,  the  coefficient  is  dependent 
on  the  mass  transfer  rate  since  high  rates  may  change  the  concentration  profiles. 
However,  for  slow  mass-transfer  rates,  such  as  those  typical  of  evaporation  processes, 
the  concentration  profile  distortion  is  negligible.  The  coeflBcient  also  depends  on  the 
position  on  the  boundary.  However,  since  the  usual  value  measured  is  the  flux  av¬ 
eraged  over  the  entire  surface  area,  this  dependence  is  neglected  so  that  measured 
values  may  be  used. 

Since  k  is  &  function  of  many  properties  of  the  fluid  A  and  gas  B,  the  ex¬ 
pression  for  a  particular  geometry  is  usually  found  in  terms  of  a  smaller  number  of 
dimensionless  parameters.  For  convective  mass  transfer  between  an  object  with  a 
characteristic  length  scale  L  and  a  surrounding  fluid,  k  is  defined  in  terms  of  the 
Sherwood  number,  ATg,  which  is  also  called  the  Nusselt  number  for  mass  transfer. 
The  Sherwood  number  is  arrived  at  through  dimensional  analysis  as  in  Bird  (1)  and 
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Clark  (2).  Previous  researchers  have  noticed  similarities  between  mass  flux  expres¬ 
sions  derived  for  different  problems  with  different  geometries: 


•  The  flux  is  proportional  to  a  characteristic  concentration  difference. 

•  The  flux  is  inversely  proportional  to  the  characteristic  length  scale,  L. 

•  The  flux  is  proportional  to  0.5  <  p  <  1. 

These  three  observations  suggest  that  multiplying  both  sides  of  (2.9)  by  Lj (pVAYa) 
yields  a  dimensionless  parameter  group.  This  group  is  the  Sherwood  number  or  the 
Nusselt  number  for  mass  transfer,  Ns'- 


^  JaqL  ^  kL 
^  pVAYa  pV 


(2.11) 


The  Sherwood  number  is  a  dimensionless  parameter  that  is  obtained  from  exper¬ 
imental  measurements.  The  form  of  the  number  is  as  a  function  of  the  Reynolds 
number.  Nr,  and  the  Schmidt  number,  Nc,  and  the  geometry  of  the  system: 


Ns  =  f{NR,  Nc,  geometry) 
where  the  Reynolds  number  is  deflned  as 

P 

and  the  Schmidt  number  is  the  deflned  as 


where  u  is  the  object’s  mean  velocity  relative  to  the  surrounding  fluid,  p  is  the  density 
of  the  fluid,  and  p  is  the  absolute  viscosity  of  the  fluid.  Notice  that  the  Schmidt 
number  has  pV  in  the  denominator.  This  is  where  the  exponent  p  comes  from  in  the 
observation  that  the  flux  is  proportional  to  pP**. 
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By  multiplying  the  Reynolds  number  by  u/L"^,  we  realize  its  physical  signifi¬ 


cance: 


Nr  = 


pu^/L  Lup 


(2.12) 


pu/LP' 

This  is  the  ratio  of  inertial  to  viscous  forces. 

We  obtain  the  mass  transfer  coefficient  by  solving  (2.11)  for  k  which  yields: 


_  p'DNs 
~  L 


(2.13) 


Now  that  we  have  the  basic  equations  needed  to  describe  evaporation,  consider 
the  evaporation  of  a  spherical  droplet  of  A  in  a  motionless  surrounding  fluid  B. 
Notice  that  in  spherical  coordinates,  the  net  bulk  flow  induced  by  ordinary  diffusion 
in  (2.10)  is  zero.  This  results  in  the  equality  of  and  jaq-  Thus: 


3aq  =  riAQ  =  k/SYA  (2.14) 

For  sufficiently  small  At,  the  mass  of  the  droplet  lost  during  time  At  is  the 
rate  of  the  mass  crossing  the  surface  times  At: 

m(t  -I-  At)  —  m(t)  =  TrD^jAo^t  (2.15) 

where  D  is  the  diameter  of  the  sphere  and  'kD^  is  the  surface  area  of  the  sphere. 

Dividing  both  side  of  (2.15)  by  At  and  taking  the  limit  as  At  0  results  in  the 

instantaneous  evaporation  rate: 

^  (2.16) 
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where  jAo  is  the  rate  of  mass  flow  per  unit  area  of  A  across  the  surface  of  the  sphere. 
Replacing  jao  with  the  expression  in  (2.14)  results  in; 


dm 

dt 


wD^kAYA 


(2.17) 


The  concentration  gradient  is  taken  to  be  AY  a  =  Yao  —  Yaooi  the  difference  between 
the  gas  phase  concentration  of  A  at  the  surface  and  the  concentration  of  A  in  the 
approaching  gas  B.  Assume  that  the  gas  concentration  of  A  at  the  surface  of  the 
droplet  is  equal  to  the  equilibrium  concentration  and  that  the  ambient  concentration 
of  A  is  zero  which  results  in  AYa  =  Yao-  Then,  substituting  k  from  equation  (2.13), 
the  change  in  mass  becomes: 


dm 

dt 


'^D^Yao 


pVNs 

D 


(2.18) 


where  the  characteristic  length  scale  is  D  and  the  empirical  relation  for  the  Sherwood 
number  for  a  sphere  surrounded  by  a  flowing  fluid  is  given  by: 


Ns  =  2.0  +  O.mNU'^Ny^ 


(2.19) 


Notice  that  the  product  p  Yao  is  the  mass  concentration  of  A  at  the  phase 
boundary,  paoj  pao  =  c  wa  Xao  where  c  is  the  total  molar  concentration,  wa  is 
the  molecular  weight  of  species  A,  and  Xao  is  the  gas  phase  mole  fraction  of  species 
A  at  the  surface.  These  quantities  are  more  easily  calculated  than  pAo-  Then,  the 
evaporation  rate  becomes: 

^  =  ttDcwaXao'DNs  (2.20) 

dt 

Prom  the  ideal  gas  law,  we  have 

c^n/V  =  P/{RT),  (2.21) 
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where  P  is  pressure,  V  is  volume,  n  is  the  number  of  moles  of  gas,  R  is  the  universal 
gas  constant,  and  T  is  the  temperature.  If  we  assume  that  A  and  B  form  an  ideal 
gas  mixture,  then  c  =  njV  —  P/  (RT)  which  yields: 

^  =  wDwaXaoVNs^  (2.22) 


If  the  gas  phase  concentration  at  the  surface  is  the  equilibrium  concentration, 
the  mole  fraction  is  equal  to  the  vapor  pressure  ,  P*',  divided  by  the  total  pressure. 
Substituting  the  expression  xao  =  Pl/^  (2.22),  the  instantaneous  evaporation 


rate  of  A  is  given  by: 


dt 


NsVwa  pi 
RT 


(2.23) 


Since  Ns  =  2  when  the  relative  droplet  velocity  is  zero,  this  expression  is  equivalent 
to  Maxwell’s  equation  for  the  evaporation  of  droplets  with  zero  velocity  with  respect 
to  the  surrounding  fluid  (see  equation  (1.14)  in  reference(5)). 


2.2.2  Heat  'Transfer.  There  are  several  mechanisms  of  heat  transport.  For 
the  evaporation  problem,  we  must  consider  heat  transfer  between  an  object  and  a 
surrounding  fluid  as  a  result  of  conduction,  convection,  and  radiation.  The  principal 
references  for  the  equations  describing  these  mechanisms  of  heat  transport  are  Bird, 
Stewart,  and  Lightfoot  (1)  and  Incropera  and  DeWitt  (7). 

When  evaporation  occurs,  heat  is  lost  through  evaporative  cooling.  Although 
the  heat  is  not  directly  transferred  to  the  surrounding  fluid,  it  is  presented  here  since 
it  effects  the  temperature  of  the  evaporating  liquid. 


2.2.2. 1  Heat  Conduction.  Heat  conduction  is  the  transfer  of  heat 
through  a  stationary  medium  due  to  a  temperature  gradient.  The  mechanism  of 
conduction  is  molecular  or  atomic  activity  where  heat  is  transferred  from  particles 
with  high  energy  to  particles  with  lower  energy. 
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Fourier’s  Law  quantifies  the  rate  of  heat  transfer  by  conduction: 


q  =  -kVT  (2.24) 

where  q  is  the  heat  flux,  k  is  the  thermal  conductivity,  and  VT  is  the  temperature 
gradient.  It  states  that  the  heat  transfer  per  unit  time  per  unit  area  is  proportional 
to  the  temperature  gradient.  Note  that  the  heat  flows  in  the  direction  of  steepest 
temperature  descent. 

Note  the  similarity  of  equation  (2.24)  with  equation  (2.7).  We  will  see  that 
the  similarities  between  mass  and  heat  transfer  will  carry  through  the  discussion 
heat  convection  as  well.  The  empirical  relationship  for  the  heat  transfer  coefficient 
is  analogous  to  that  for  the  mass  transfer  coefficient  where  the  Sherwood  number  is 
replaced  by  the  Prandtl  number  and  the  pV  product  replaced  by  the  thermal  conduc¬ 
tivity.  In  fact,  the  mass  transfer  coefficient  is  often  deduced  by  making  an  analogy 
to  the  corresponding  heat  transfer  coefficient  and  performing  the  replacements  just 
stated. 


2. 2. 2. 2  Heat  Convection.  Heat  convection  is  the  heat  transfer  be¬ 
tween  an  object  and  a  flowing  medium.  It  is  composed  of  two  mechanisms,  conduc¬ 
tion  and  advection.  For  example,  consider  a  fluid  flowing  over  a  heated  surface  as  in 
figure  2.2.  Because  of  friction  between  the  fluid  and  the  surface,  the  velocity  of  the 
fluid  varies  from  0  to  « throughout  a  region  called  the  hydrodynamic  boundary  layer. 
The  temperature  also  varies  in  a  region  called  the  thermal  boundary  layer,  from  Tg 
at  the  surface  to  Too,  the  ambient  fluid  temperature.  The  heat  is  moving  away  from 
the  surface  because  the  surface  is  being  heated  to  a  higher  temperature  than  the 
fluid.  It  is  a  result  of  both  conduction  to  the  boundary  layer  due  to  a  temperature 
gradient,  and  the  bulk  flow  carrying  heated  fluid  away. 
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Flow 


Velocity,  u(y)  Temperature,  T(y) 


Figure  2.2  Velocity  and  temperature  profiles  in  the  heat  transfer  to  a  fluid  flowing 
over  a  heated  surface. 
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There  are  two  limiting  types  of  heat  convection,  free  and  forced  convection. 
Although  most  processes  are  placed  into  one  of  the  two  categories,  it  is  possible  for 
both  types  of  convection  to  be  significant  simultaneously. 

Forced  convection  occurs  when  the  object  is  in  motion  with  respect  to  the 
surrounding  medium,  as  in  the  example  of  figure  2.2.  The  fiuid  is  heated  by  the 
object  and  then  carried  away  by  some  external  force  such  as  gravity  pulling  an 
object  through  the  air  or  a  fan. 

Free  convection  occurs  when  the  surrounding  fluid  is  set  into  motion  by  density 
differences.  As  the  fluid  is  heated  by  the  object,  it  becomes  less  dense.  The  medium 
will  then  move  due  to  increased  buoyancy.  The  result  is  a  net  transport  of  heat  away 
from  the  object  called  free  convection.  In  the  example  of  figure  2.2,  if  u  — >•  0,  then  we 
would  have  free  convection.  The  density  of  the  fluid  above  the  surface  would  begin 
decrease  as  it  was  heated  by  the  surface.  This  would  cause  the  fluid  near  the  surface 
to  rise  and  be  replaced  with  cooler  fluid.  As  the  heated  fluid  rises,  it  transports  heat 
away  from  the  surface. 

For  either  type  of  convection,  Newton’s  law  of  cooling  relates  the  convective 
heat  flux  to  the  temperature  gradient; 

Qe  =  hVT  (2.25) 

where  qc  is  the  convective  heat  flux,  VT  is  the  temperature  gradient,  and  h  is  the 
convection  heat  transfer  coefficient.  The  heat  transfer  coefficient  depends  on  the 
boundary  layer  conditions,  which  are  influenced  by  the  object’s  geometry  and  sev¬ 
eral  of  the  fluid  properties  such  as  the  viscosity,  density,  heat  capacity,  etc..  For 
complex  situations  involving  object-fluid  boundaries  and  convective  transport,  the 
heat  transfer  coefficient,  like  the  mass  transport  coefficient,  must  be  fitted  to  exper¬ 
imental  data.  By  making  an  analogy  to  convective  mass  transport,  the  empirical 


2-12 


relationship  for  the  heat  transfer  coefficient  is  defined  as: 

h  =  (2.26) 

where  Nu  is  the  Nusselt  number  for  heat  transfer,  and  L  is  the  characteristic  length 
of  the  object.  The  empirical  relationship  for  Nu  is  analogous  to  that  for  Ns'. 

Nu  =  /{Nr,  Np,  geometry) 

where  Np  =  fiCp/n  is  the  Prandtl  number,  the  heat  transfer  analog  to  the  Schmidt 
number,  fj,  is  the  fiuid  viscosity  at  the  surface  and  Cp  is  the  gas  phase,  constant 
pressure  heat  capacity  of  the  fiuid  at  the  surface. 

For  convective  heat  flux  between  a  spherical  droplet  and  a  surrounding  fluid, 
the  Nusselt  number  is  given  by  a  relationship  similar  to  that  for  the  Schmidt  number: 

Nu  =  2.0  +  0.6N]I^nI/^ 

If  we  also  assume  that  we  have  spherical  symmetry,  then  VT  =  dT/ dr  which  is  set 
equal  to  T*  -  Too.  The  heat  flux  normal  to  the  droplet  surface  from  the  surrounding 
fluid  is  then: 

*  =  5(2.0  +  OSN'J^N'J^XT,  -  T^)  (2.27) 

where  Qc  is  the  radial  heat  flux  through  conduction  with  convection  and  the  droplet 
diameter  D  is  the  characteristic  length  scale. 

2. 2. 2. 3  Heat  Radiation.  For  the  range  of  magnitudes  for  the  temper¬ 
ature  gradient  between  the  atmosphere  and  either  a  released  droplet  or  the  resultant 
film  of  fuel,  heat  transfer  by  radiation  will  be  small  relative  to  convective  trans¬ 
fer.  However,  it  is  presented  here  for  completeness  because  it  is  incorporated  into 
previous  work. 
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The  principles  of  radiative  transfer  are  derived  from  the  analysis  of  black  bod¬ 
ies,  objects  for  which  absorptivity  is  unity  over  all  wavelengths  and  temperatures. 
The  absorptivity  a  is  defined  as  the  ratio  5“/ g*,  where  g®  is  the  fiux  of  energy  ab¬ 
sorbed  and  g*  is  the  incident  fiux  of  energy.  Any  real  body  will  have  absorptivity 
less  than  unity. 

The  total  energy  emitted  by  a  black  body,  gf,  is  given  by  the  Stefan-Boltzmann 

law: 

=  aT^  (2.28) 

where  a  is  the  empirically  derived  Stefan-Boltzmann  constant  and  T  is  the  tempera¬ 
ture  in  Kelvin.  With  the  exception  of  fiuorescent  bodies,  the  energy  emitted  by  real, 
nonblack  bodies  will  be  less  than  indicated  in  (2.28).  How  much  less  is  characterized 
by  the  emissivity  of  the  object  e,  which  is  defined  as: 

e  =  ^  <  1  (2.29) 

It 

where  g®  without  the  subscript  b  is  the  energy  emitted  by  the  nonblack  body.  There¬ 
fore,  the  energy  emitted  by  a  real  body  is  given  by: 

g®  =  eaT^  (2.30) 

To  quantify  the  radiant  energy  transfer  between  two  bodies,  we  must  consider 
only  the  radiation  directed  at  the  bodies.  To  see  this  is  not  necessarily  the  same  as 
the  total  energy  emitted,  consider  the  transfer  between  two  spheres  separated  by  a 
vacuum.  Although  energy  is  radiated  from  the  entire  surface  of  a  body,  only  that 
which  is  directed  at  the  other  body  contributes  to  any  exchange. 

When  an  object  is  surrounded  by  a  fluid,  however,  all  of  the  energy  radiating 
from  the  object  is  directed  toward  the  surrounding  fluid  and  we  can  simply  use 
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(2.30).  The  net  flux  of  energy  from  the  object  to  the  surrounding  fluid  is  then: 

5,  =  -  a,Ti,)  (2.31) 

where  e*  is  the  emissivity  of  the  object,  Tg  is  the  temperature  of  the  object,  is  the 
absorptivity  of  the  object,  and  Too  is  the  temperature  of  the  fluid. 

2.2.24  Evaporative  Cooling.  There  is  a  certain  amount  of  heat  en¬ 
ergy,  called  the  heat  of  formation,  associated  with  each  of  the  three  phases,  solid, 
liquid,  and  gas.  When  ‘a  substance  moves  from  one  phase  to  another,  there  is  a 
change  in  heat.  The  change  in  heat  that  occurs  when  a  substance  moves  from  liquid 
to  gas  is  quantifled  by  the  latent  heat  of  vaporization,  AHyap-  The  latent  heat  of 
vaporization  gives  the  amount  of  heat  lost  per  unit  mass  of  the  substance  that  is 
vaporized.  Thus,  the  flux  of  heat  due  to  evaporative  cooling,  Qh  is  given  in  term  of 
the  change  in  mass  as: 

AHyg^p  dm  fiy 

*  =  —n 

where  A  is  the  surface  area  of  the  evaporating  liquid  which  is  in  the  denominator 
in  order  to  be  consistent  with  units  of  the  previous  modes  of  heat  flux,  conduction, 
convection,  and  radiation. 
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2.3  Existing  Droplet  Evaporation  Models 

2.3.1  Overview.  The  evaporation  of  jettisoned  jet  fuels  has  been  studied 
for  several  decades.  In  particular,  two  efforts  have  led  to  the  publishing  of  fuel 
evaporation  descriptions  to  predict  the  fate  of  jettisoned  fuel  from  aircraft.  Each 
investigation  is  based  upon  different  approaches  to  the  problem. 

The  first  approach  was  developed  by  Lowell  at  NASA(IO).  It  can  be  shown 
that  Lowell’s  approach  is  based  on  the  heat  and  mass  transfer  principles  presented 
in  section  2.2.  Lowell’s  work  was  further  improved  by  Clewell(3),  followed  by  Pfeif¬ 
fer,  Quinn,  and  Dunge  at  AFIT(13).  The  advances  in  the  dispersion  modeling  of 
jettisoned  fuel  are  presented  here  because  they  must  be  used  to  determine  the  initial 
amount  of  fuel  reaching  the  ground  crew. 

The  second  approach  of  Runge  et  al.  (16)  is  based  on  Williams’ (18)  and 
Law’s(9)  work  on  droplet  combustion. 

The  assumptions  and  equations  of  the  two  approaches  are  presented  in  this 
chapter.  As  they  are  presented,  the  equations  of  section  2.2  are  referenced  to  demon¬ 
strate  the  scientific  foundations  of  the  models.  This  knowledge  of  the  principles  in¬ 
corporated  will  make  it  possible  to  determine  what  modifications  must  be  made  to 
extend  the  equations  to  describe  the  evaporation  of  fuel  from  a  non-spherical  surface. 

2.3.2  Lowell,  Clewell,  and  Pfeiffer.  Lowell  developed  a  description  of 
the  evaporation  and  dispersion  of  jettisoned  JP-4  (10)  and  JP-1  (11)  jet  fuel.  For 
the  JP-4  model,  the  droplet  temperature  was  assumed  to  always  equal  the  ambient 
temperature  to  simplify  the  calculation.  For  his  later  research  on  the  less  volatile 
JP-1,  he  included  a  droplet  cooling  calculation  and  an  adaptive  time  step.  The 
adaptive  time  stepping  scheme  increased  the  time  step  if  the  evaporation  rate  was 
small.  Lowell  established  upper  bounds  on  the  amount  of  mass  evaporated  and 
distance  fallen  during  a  time  step.  If  the  upper  bounds  are  exceeded,  the  time  step 
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is  shortened  to  maintain  the  accuracy  of  the  approximate  solution.  This  scheme 
significantly  reduces  the  computational  time  of  the  simulation. 

Jet  fuel  is  a  very  complex  mixture  of  hydrocarbons  and  to  characterize  the 
evaporating  substance  perfectly  is  not  practical.  Furthermore,  variations  in  the 
refining  process  result  in  variations  in  the  composition  of  the  fuel.  Therefore,  the  fuel 
has  been  characterized  by  a  mixture  of  a  finite  number  of  species  that  approximate 
the  physical  behavior  of  the  actual  compounds  in  the  mixture.  The  individual  fuel 
species  are  characterized  by  the  volume  fraction,  molecular  weight,  boiling  point,  and 
reference  density  at  20  degrees  Celsius.  The  characterization  developed  by  Lowell 
used  only  a  10  component  characterization  of  the  fuels. 

Lowell  also  developed  a  dispersion  model  to  determine  the  concentration  of  the 
fuel.  He  characterized  the  jettisoned  plume  as  an  infinite  line  source  and  calculated 
a  TniniTTinTn  altitude  for  jettisoning  at  which  the  fuel  presented  no  fire  hazard.  Lowell 
noted,  however,  that  further  study  was  needed  to  predict  ground  contamination  by 
the  jettisoned  fuel. 

Clewell’s  study  (3)  is  an  improvement  of  the  work  of  Lowell  to  predict  the  evap¬ 
oration  of  jettisoned  aircraft  fuel.  Clewell’s  model  is  refined  by  the  inclusion  of  ex¬ 
perimental  data  provided  by  the  Arnold  Engineering  Development  Center  (AEDC), 
an  initial  droplet  temperature  estimate  based  on  the  speed  of  the  jettisoning  air¬ 
craft,  and  a  more  detailed  fuel  characterization  with  33  components  for  JP-4  and  27 
components  for  JP-8.  The  more  complete  fuel  characterization  allows  the  model  to 
simulate  99.9%  evaporation  of  the  initial  mass. 

Clewell  developed  a  simpler  dispersion  model  than  Lowell’s  to  use  with  the 
evaporation  model  to  predict  ground  contamination  by  the  fuel.  Although  his  dis¬ 
persion  box  model  was  less  accurate  than  Lowell’s  infinite  line  source,  Clewell  was 
able  to  use  it  to  get  an  upper  bound  on  the  amount  of  ground  contamination. 
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Pfeiffer’s  work,  while  at  the  Air  Force  Institute  of  Technology  (AFIT),  was 
based  on  the  evaporation  equations  of  Lowell  and  Clewell.  He  improved  the  de¬ 
scription  by  with  a  more  accurate  Gaussian  dispersion  model  and  by  incorporating 
of  actual  upper  atmosphere  data  obtained  from  several  hundred  stations  worldwide 
that  report  vertical  profiles  of  temperature,  humidity,  pressure,  and  wind  velocity. 
This  more  accurate  atmosphere  characterization  leads  to  major  improvement  in  the 
model’s  accuracy  as  far  as  fuel  jettisoning  is  concerned.  This  more  accurate  disper¬ 
sion  model  is  a  significant  improvement  in  predicting  the  concentration  of  fuel  upon 
reaching  the  ground  crew  (13,  12). 

2. 3. 2.1  Assumptions. 


1. 

2. 

3. 

4. 

5. 

6. 

7. 


Quasi-steady  state  evaporation.  The  droplet  diameter  is  fixed  at  each  time 
step  and  then  updated  after  the  new  mass  and  temperature  are  calculated. 

The  droplet  falls  at  the  terminal  velocity  for  its  current  diameter,  density,  and 
altitude. 

Raoult’s  law  applies. 

Each  component  evaporates  independently  of  the  rest  of  the  mixture. 

The  initial  fuel  temperature  is  equilibrated  with  the  temperature  of  the  skin 
of  the  aircraft. 

The  atmosphere  is  standard.  This  assumption  does  not  have  to  be  made  if  the 
actual  profiles  are  obtained. 

The  mixture  is  ideal  so  that  the  droplets  volume  is  the  sum  of  the  volumes 
of  the  components.  This  assumption  closely  approximates  the  actual  solution 
because  all  of  the  components  of  the  mixture  are  hydrocarbons. 
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2. 3. 2. 2  Model  Equations.  Following  Clewell  (3) ,  the  evaporation  rate 
of  the  component  is: 


drui 

dt 


^D%py, 


i  =  l..N 


(2.33) 


where  ki  is  the  mass  transfer  coefficient  of  component  i  and  is  given  by: 


DRT^ 


(2.34) 


Notice  that  equation  (2.33)  is  equivalent  to  equation  (2.23)  when  we  subscript  the 
properties  for  a  multicomponent  mixture  and  thereby  let  the  diffusivity  “Di  denote 
the  multicomponent  diffusivity. 

When  we  assume  Raoult’s  Law  applies,  the  vapor  pressure  exerted  by  compo¬ 
nent  z  in  a  mixture  is  calculated  using: 


py = xi,epr 


(2.35) 


where  is  the  pure  liquid  vapor  pressure  of  component  i,  Py  is  the  pressure 
exerted  by  component  z  in  a  mixture,  and  Xi,e  is  the  mole  fraction  of  species  z  in  the 
liquid  phase. 


In  equation  (2.34),  the  subscript  on  the  Sherwood  number  is  necessary  because 
it  is  a  function  of  the  Schmidt  number,  which  is  a  function  of  the  diffusivity  of  the 


z*'*  component,  Vi". 


N,a  = 


JL 

pVi 


(2.36) 


The  diffusivity  of  component  z  in  air  is  computed  following  Clewell  who  derived  the 
expression  from  one  given  in  Bird  et  al.  (1): 


2.66  *  io-^r^(iM  + 

^.^(^‘'’+0.31)= 


(2.37) 
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where  Tm  is  the  mean  of  the  droplet  and  ambient  temperatures,  14  is  the  molar 
volume  of  component  i  at  its  normal  boiling  point,  Patm  and  the  subscript  a  denotes 
air. 

The  differential  equation  is  solved  using  a  finite  difference  approximation  of 
the  time  derivative  to  obtain  a  mass  for  the  component  of  the  fuel: 

Ami  =  T^D^kiPiXi/^t,  i  =  (2.38) 

where  Arui  =  mi{t  +  At)  —  mi{t),  and  A<  is  the  time  step  for  the  approximation 
which  is  adjusted  as  in  Lowell’s  work(ll). 

The  droplet’s  net  rate  of  change  in  thermal  energy,  Q,  is  computed  by  per¬ 
forming  an  energy  balance: 

Q  =  Heat  in  per  unit  time  -  Heat  out  per  unit  time 

Recall  from  equation  (2.31)  in  section  2.2  that  the  fiux  of  heat  transferred  to 
the  droplet  by  radiation  from  the  atmosphere  is  given  by: 


Qr  =  a{adT^  -  erfTj) 

where  the  subscript  d  denotes  a  property  of  the  droplet. 

Prom  equation  (2.27)  in  section  2.2,  the  flux  of  heat  transferred  to  the  droplet 
by  conduction-convection  is: 


qc  =  h  AT  =  h  (Too  -  Td) 


The  Prandtl  number  and  thermal  conductivity  used  in  the  calculation  of  the  heat 
transfer  coefficient  are  assumed  to  be  approximately  those  of  air  (3). 
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In  addition,  there  will  be  a  flux  of  heat  loss  by  the  droplet  due  to  evaporative 
cooling  that  is  given  by  equation  (2.32)  in  section  2.2: 

A  W 

L^riyap  at 


where  dm/dt  is  the  sum  of  the  species  evaporation  rates,  drrii/dt. 

Furthermore,  if  the  fuel  is  jettisoned  at  high  altitude,  the  droplet  will  be  heated 
signiflcantly  by  solar  radiation.  The  flux  of  heat  to  the  droplet  by  this  mechanism 
is  given  by: 

L{1  -  a) 

4 

where  a  is  the  droplet’s  albedo  and  L  is  the  solar  insolation  rate  at  altitudes  typical 
of  fuel  jettisoning  events.  Clewell  uses  a  rather  high  value  of  1000  for  the  solar 
insolation  rate  to  account  for  the  high  altitudes  encountered.  This  value  is  not 
adjusted  as  the  droplet  falls. 

Combining  these  relations  into  an  energy  balance  by  adding  the  heat  trans¬ 
ferred  to  the  droplet  and  subtracting  the  heat  loss  by  evaporative  cooling  results  in 
the  rate  of  heat  transferred  to  the  droplet: 

Q  =  TrD^{qr  +  Qs  +  Qc  ~  Qh)  (2.39) 

where  is  the  surface  area  of  the  droplet  and  qh  is  subtracted  because  it  has  been 
formulated  in  terms  of  the  positive  evaporation  rate  instead  of  the  negative  change 
in  mass.  Since  Q  is  the  change  in  energy  per  unit  time,  Q  divided  by  m  ci  results  in 
the  change  in  temperature  with  respect  to  time: 


dTj  _  Q 

dt  met 


(2.40) 
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The  droplet  surface  temperature,  Ta,  is  then  calculated  in  one  of  two  ways 
depending  on  whether  the  droplet  is  warmer  or  cooler  than  the  ambient  tempera¬ 
ture.  If  the  droplet  is  warmer  than  the  ambient  temperature,  the  droplet  cooling  is 
calculated  by  using  a  forward  difference  approximation  to  the  derivative: 


ATrf  = 


dt 


At 


If  the  droplet  is  cooler  than  the  ambient  air,  a  steady-state  calculation  is  performed. 
The  difference  between  the  ambient  and  droplet  temperature  is  calculated  by  setting 
Q  =  0  and  solving  for  the  equilibrium  temperature  difference  at  the  conditions  for 
the  current  time  step.  Since  this  equation  is  nonlinear,  a  Newton-Raphson  method 
is  used  to  solve  for  the  temperature  iteratively: 


Tn+l  Tn 

where  the  prime  denotes  differentiation  with  respect  to  T.  The  initial  guess  is  To  = 
T{t)  and  the  iteration  continues  until  T^+i  -  Tn  <  0.01.  If  the  temperature  change 
calculated  from  (2.41)  is  too  great,  (2.38)  is  recalculated  using  the  new  temperature. 
This  process  is  repeated  until  a  steady  state  temperature  is  reached. 

The  droplet’s  altitude  is  updated  by 


H{t  -I-  At)  =  H{t)  -  AH  =  Hit)  -  uit)At  (2.42) 

where  tt(t)  is  the  droplet’s  terminal  velocity  for  its  size  at  time  t. 

This  algorithm  incorporates  an  adaptive  time  step.  The  time  step  is  adjusted 
based  on  the  magnitude  of  the  changes  in  mass  and  altitude.  If  either  of  the  changes 
calculated  in  equations  (2.38)  and  (2.42)  are  greater  than  their  prescribed  maximums, 
At  is  halved  and  Am,  AH,  and  AT  are  recalculated  until  the  changes  in  the  mass 
and  altitude  are  below  their  maximums. 
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Once  AH  and  Am  are  within  acceptable  limits,  droplet  conditions  are  updated 
for  use  in  the  next  iteration.  If  AH  and  Am  are  now  below  their  prescribed  mini- 
mums,  the  time  step  is  doubled  and  the  calculation  continues  using  the  new  droplet 
and  ambient  conditions  as  input  for  the  next  time  step. 

The  droplet  diameter  is  calculated  from  the  new  droplet  volume  as  follows: 

Vi  =  m/PiKT)  (2.43) 

Vi  =  Y,V  (2.44) 

i=l 

D  =  (2.45) 

TT 

where  p{T)  is  a  density  correction  factor  to  account  for  the  temperature  dependence 
of  p,  Vi  is  the  volume  that  species  i  contributes  to  the  total  volume,  V^. 

2.S.3  Runge,  Teske,  Polymeropolous.  Two  models  are  considered  in  this 
section,  both  of  which  are  based  on  Law’s  rapid  mixing  model  approach  to  describing 
combustion  of  fuels  (9)  as  discussed  in  this  section.  The  first  model  is  found  in  a 
master’s  thesis  by  Runge  (17).  The  second  was  presented  in  an  article  by  Runge, 
Teske,  and  Polymeropolous  (16). 

To  burn  a  single  droplet  of  fuel,  the  condensed  fuel  must  first  be  vaporized. 
Therefore,  the  rate  at  which  a  droplet  burns  is  controlled  by  the  rate  of  evaporation. 
For  this  reason,  combustion  models  include  a  calculation  of  rates  of  evaporation.  The 
system  described  here  assumes  a  single,  multicomponent  fuel  droplet  evaporating  in 
a  quiescent  atmosphere.  The  rate  obtained  through  this  analysis  is  then  corrected 
for  eflfects  of  convection  using  the  methods  as  discussed  in  section  2.2. 

2.3.3. 1  Assumptions. 

1.  Spherical  symmetric  droplets  with  empirically  corrected  vaporization  enhance¬ 
ment  due  to  droplet  motion. 
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2.  Spatially  uniform,  temporally  varying  temperature  and  component  concentra¬ 
tions  in  the  droplet. 

3.  Raoult’s  law  applies. 

4.  The  thermal  conductivity,  A,  and  gas  and  liquid  phase  specific  heats,  Cp  and 
Cf,  are  constant  with  respect  to  time  and  fuel  component. 

5.  The  Lewis  number,  Nl  =  NjNp,  is  unity  for  the  gas  phase  which  implies  that 
the  Schmidt  and  Prandtl  numbers  are  equal.  Here,  the  Schmidt  and  Prandtl 
numbers  are  both  set  equal  to  one.  This  assumption  also  implies  that  the 
thermal  and  mass  diffusivities  are  equal. 

6.  The  value  of  pT>im  =  pV,  the  density  times  mass  diffusivity,  is  the  same  for  all 
gaseous  species. 

7.  The  gas  and  liquid  mixtures  behave  ideally. 

Runge  states  in  his  thesis  that  assumptions  5  and  6  are  not  justified  and  are 
used  for  convenience.  Using  constant  properties  is  not  necessary  to  obtain  a  solution, 
but,  it  simplifies  the  calculations.  These  simplifications  are  taken  into  account  by 
fitting  parameters  to  experimental  data  as  described  next. 

The  more  volatile  fuel  components  will  evaporate  faster  than  the  less  volatile 
components.  As  the  mixture  evaporates,  the  less  volatile  components  increase  in 
concentration,  changing  the  fuel  properties.  These  changes  will  cause  erroneous  re¬ 
sults  from  a  calculation  of  multiple  component  evaporation  using  constant  properties 
that  are  the  same  for  all  species.  Runge  et  al.  remedy  this  problem  with  a  weighted 
empirical  correlation  of  pV  based  on  the  percent  volume  of  liquid  remaining: 


pV  =  VppVf,,  -f-  (1  -  Vp)pVi,,  (2.46) 

where  Vp  is  the  percent  volume  of  liquid  remaining,  Vhi  is  the  upper  limit  diflFusivity 
constant,  Vio  is  the  lower  limit  diffusivity  constant.  These  constants  are  determined 
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through  empirical  parameter  fitting  to  experimental  data.  Time  is  then  scaled  by 
pVfpVr,  where  the  subscript  r  denotes  a  constant  reference  value,  to  correct  for  the 
changes  in  the  actual  pT>  product. 

The  descriptions  presented  in  Runge’s  thesis  (17)  and  the  article  (16)  report 
different  values  for  the  parameters  T>hi  and  T>io-  Reference  (16)  obtains  the  best-fit 
parameters  shown  in  table  2.1.  However,  reference  (17)  reports  that  the  best  fit 
occurs  when  the  values  are  those  found  in  table  2.2. 


Fuel 

Vio  *  10^,  cm^/s 

Vhi  *  10®,  crri^/s 

JP-4 

2.0 

4.0 

JP-8 

2.0 

11.0 

Table  2.1  pT>  coefficients  reported  in  the  article. 


Fuel 

Vio  *  10®,  cm^/s 

Vhi  *  10®,  cm^/s 

JP-4 

2.0 

7.0 

JP-8 

4.0 

14.0 

Table  2.2  pV  coefficients  used  in  Runge’s  thesis. 

2. 3. 3.2  Equations.  We  begin  with  equation  (2.8)  in  section  2.2,  Pick’s 
first  law  of  diffusion: 


t^a  =  PaVa  =  YaIxia  +  ns)  -  p'Dab'^Ya 


If  we  assume  spherical  symmetry,  then  this  becomes  a  one-dimensional  problem  and 
the  flux  of  A  is  given  by: 

ua  =  YA{nA  +  n^)  -  pVab-^  (2.47) 

For  multicomponent  mixtures,  an  effective  binary  diffusivity  Vim  for  the  dif¬ 
fusion  of  i  in  a  mixture  of  N  species  can  be  defined  by  an  analogous  equation  as  in 
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Bird  et  al.  (1): 

N  fiY. 

n,  =  Yi'£’^i-p-Di„Y  (2.48) 

where  the  subscript  m  denotes  the  rest  of  the  mixture.  Multiplying  both  sides  of 
equation  (2.48)  by  47rr^  results  in  an  expression  for  the  mass  of  species  i  leaving  the 
total  area  of  a  sphere  of  radius  r,  Mf. 

N  dY- 

k=l 

Then,  dividing  both  sides  by  47rropI>im  yields  equation  (1)  of  Law: 

=  (2.49) 

where  A4=  M / {ATtr^p'Di^  is  the  non-dimensional  mass  evaporation  rate  and  R= 
r/ro. 

The  equations  for  the  changes  in  radius,  mass,  and  temperature  are  all  based 
upon  the  normalized  aggregate  vaporization  rate,  Mf  =  as  will  be  de¬ 

scribed  in  section  2. 3.3.2. 

The  energy  equation  used  by  Law  (9)  can  be  arrived  at  from  the  more  familiar 
form  of  the  energy  equation  (see  equation  (6-105)  of  Kuo  (8)): 

=  ^)  +  (2.50) 

dr  dr  Cp  dr 

where  A  is  the  thermal  conductivity,  and  Q  is  the  heat  of  reaction  per  unit  mass  of  fuel 
during  combustion.  For  pure  evaporation  (no  combustion),  Q  is  zero.  Following  the 
analysis  presented  by  Kuo,  and  then  non-dimensionalizing  the  variables,  we  arrive 
at  equation  (2)  of  Law: 

MpiT  -  Ts)  -  =  -Mf{C  +  H)  (2.51) 
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where  T  =  Tcp/Lr,  Cp  is  the  constant  gas  specific  heat,  Lr  is  a  reference  latent 
heat  of  vaporization,  and  C  =  Eili  is  the  species  weighted,  non-dimensional 

latent  heat  of  vaporization.  The  heat  input  function,  H,  is  the  difference  between  the 
heat  transferred  to  the  droplet  and  that  needed  for  vaporization.  It  is  not  explicitly 
defined  and  is  obtained,  along  with  equations  (2.52), (2.55),  and  (2.57),  by  solving 
equations  (2.49)  and  (2.51). 


The  equation  of  change  for  the  non-dimensional  mass  of  species  i  is  given  in 
reference  (17)  by: 


drui 

dt 


=  —M-iR  i  =  l..A^, 


(2.52) 


where  Mi  =  eM  is  the  spherically  symmetric  mass  vaporization  rate  and  R  =  r/ro 
is  the  non-dimensional  radius  with  r  being  the  radius  at  time  t  and  ro  the  initial 
radius  of  the  droplet.  The  spherically  symmetric  mass  vaporization  rate,  M,  is  given 
as  in  Law  (9): 

M  =  ln{l  +  B)  (2.53) 


where 

B  =  ^  (2.54) 

is  the  driving  force  for  evaporation,  the  mass  transfer  or  Spalding  number.  In  calcu¬ 
lating  Mi,  Ci  =  Yi/Yp  is  the  ratio  of  the  mass  fraction  of  i  in  the  vapor  phase  at  the 
surface  of  the  droplet,  taking  into  account  the  presence  of  air,  to  the  sum  of  those 
mass  fractions,  Yp. 


The  vaporization  rate  is  empirically  adjusted  for  convective  effects  using  the 
Sherwood  number,  which  is  computed  using  equation  (2.19)  of  section  2.2: 


Ns  =  2.0  + 


where  the  Schmidt  number  has  been  set  equal  to  one.  Notice  that  since  it  has  been 
assumed  that  Nc  =  Np,  the  Sherwood  number  is  the  same  as  the  Nusselt  number 
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for  heat  transfer.  As  we  will  see  next,  M  is  also  used  in  the  calculation  of  the  heat 
transfer. 

The  equation  for  the  change  in  temperature  with  respect  to  time  is: 


dT  HM 
dt  p(,B? 


(2.55) 


where  H  is  the  heat  input  function,  M  is  the  spherically  symmetry  mass  vaporization 
rate  of  the  droplet,  and  pi  is  the  density  of  the  liquid  fuel. 

The  heat  input  function  is  calculated  by 


H  =  {I-  Yf){Too  -  Ts)IYf  -  A  (2.56) 


where  Too  =  TooCplLr  is  the  non-dimensional  ambient  temperature,  Ts  =  T^Cp/Lr 
is  the  non-dimensional  droplet  temperature,  Cp  is  the  gas  phase  specific  heat  at  the 
surface  of  the  droplet,  C  =  E"=i  A/A  is  the  latent  heat  of  vaporization,  Li  is  the 
latent  heat  of  species  i,  and  Lr  is  a  reference  latent  heat. 

The  equation  for  the  change  in  non-dimensional  radius  is  given  by: 


m  _  -M 
dt  ZRq 


(2.57) 


where  q  =  pi/pr  is  the  liquid  density  divided  by  a  reference  density.  If  this  equation 
is  multiplied  by  ZB?q,  then  we  obtain: 


=  ^qR?  =  -MR 
dt  dt 

since  g  is  assumed  to  be  constant.  This  final  form  is  the  equation  implemented  in 
Range’s  code  that  is  listed  in  the  appendix  of  his  thesis.  Runge  then  updates  R^g 
and  recovers  R  from  this  value. 


2-28 


The  equations  presented  in  the  article  by  Runge  et.  al  (16),  where  most  of  the 
experimental  data  in  chapter  III  appears,  are  different  from  those  found  in  Runge’s 
thesis  (17).  The  non-dimensional  spherically  symmetric  mass  vaporization  rate  is 
given  in  reference  (16)  by: 

M  =  RpVln(l  -I-  B)  (2.58) 

This  rate  is  then  corrected  to  account  for  convective  effects  by  multiplying  by  the 
Sherwood  number,  which  is  given  by: 

Ns  =  (2.0  -f  0.6iVi/2Ar]/3)(i  + 

In  the  case  of  heat  transfer,  it  is  corrected  for  convective  effects  by  multiplying  by 
the  Nusselt  number,  which  is  given  by: 

K  =  (2  +  0.6Ari/"iVy")(l  -h  Too  - 


The  time  rate  of  change  of  non-dimensional  radius,  R,  temperature,  %,  and  mass, 
m,  is  then  given  by: 


dt 

drfii 

dt 


M 

ZRq 

-CiM 


dTs  _  CpHM 
dt  ctR?Q 


(2.59) 

(2.60) 
(2.61) 


where  q  is  the  liquid  specific  heat  of  the  mixture.  The  heat  input  function,  H  is 
given  by: 

N  N 

H  =  Q'f-C'f  (2.62) 

where  C  is  given  as  in  reference  (17),  and  Q  is  given  by: 


Q  = 


'Too  T s 

B 
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Unfortunately,  we  were  not  able  to  use  the  equations  presented  in  reference  (16) 
to  reproduce  exactly  the  predictions  presented  in  that  paper.  This  will  be  further 
discussed  in  chapter  III. 
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III.  Comparison  of  the  Droplet  Evaporation  Calculations 

5.1  Overview 

In  this  chapter,  evaporation  calculations  of  Clewell  (3),  Runge  et  al.  (16),  and 
Runge’s  thesis  (17)  are  compared.  Clewell’s  algorithm  was  chosen  because  it  is  more 
advanced  than  Lowell’s  and  the  contribution  by  Pfeiffer  et  al.  was  improvement 
to  the  dispersion  calculations  which  are  not  considered  here.  The  evaporation  al¬ 
gorithms  of  Clewell  and  Pfeiffer  should  predict  almost  identical  results  and  a  short 
comparison  of  the  predictions  is  performed  to  verify  that  this  is  the  case. 

5.2  Comparison  of  the  Predictions  of  Clewell  and  Pfeiffer 

Pfeiffer’s  algorithm  was  available  in  electronic  format  as  part  of  his  thesis  at 
AFIT  (12).  Clewell’s  algorithm  was  implemented  from  the  code  listing  found  in 
reference  (3). 

The  models  of  Clewell  and  Pfeiffer  should  be  almost  identical  with  respect 
to  evaporation  because  Pfeiffer  used  the  model  of  Clewell  to  predict  evaporation 
rates,  focusing  instead  on  improving  the  dispersion  model.  Figure  3.1  compares  the 
predictions  of  the  evaporation  of  JP-4.  The  initial  temperature  was  calculated  from 
an  aircraft  velocity  of  175  m/s,  the  ground  temperature  was  273.15  K,  the  initial 
diameter  was  500  pm,  and  the  jettison  altitude  was  1500.0  m.  Notice  from  figure 
3.1  that  there  is  difference  in  the  diameter  predictions  of  Clewell  and  Pfeiffer. 

The  major  reason  is  the  difference  in  the  temperatures  predicted.  Clewell’s 
algorithm  predicts  a  large  decrease  in  temperature  initially.  This  is  due  to  evapo¬ 
rative  cooling  and  Pfeiffer’s  predictions  should  be,  at  least  qualitatively,  exhibiting 
this  same  behavior.  Although  the  two  algorithms  are  supposed  to  be  the  same,  the 
temperature  predicted  by  Pfeiffer  does  not  match  the  prediction  of  Clewell’s  algo¬ 
rithm,  as  shown  in  figure  3.2.  This  would  explain  why  Clewell’s  algorithm  predicts 
a  different  diameter  than  Pfeiffer’s  algorithm.  Pfeiffer’s  temperature  prediction  is 
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Figure  3.1  Comparison  of  Clewell’s  and  Pfeiffer’s  predictions  of  JP-4  droplet  diam¬ 

eter  versus  time  during  evaporation.  Jettison  velocity  =  175m/s,  Too  = 
0"  C  at  the  ground,  Dq  =  500  /xm,  jettison  altitude  =  1500  m. 


initially  higher  than  Clewell’s,  leading  to  more  evaporation  and  a  smaller  droplet 
diameter  as  shown  in  figure  3.1. 

In  the  version  of  Pfeiffer’s  algorithm  that  was  provided  with  the  AFIT  thesis, 
the  calculation  of  the  steady  state  temperature  is  erroneous.  The  calculation  of  Q 
and  Q '  for  the  Newton-Raphson  iterative  method  had  to  be  corrected.  Furthermore, 
the  droplet  temperature  was  not  being  updated  after  the  steady-state  temperature 
was  calculated  as  indicated  by  the  short  intervals  where  the  predicted  temperature 
does  not  change  in  figure  3.2.  Thirdly,  the  algorithm  only  computed  the  droplet 
cooling  at  the  first  time  step  and  then  switched  to  strictly  computing  the  steady 
state  temperature.  The  initial  calculation  of  the  number  of  moles  of  fuel  also  had 
to  be  corrected.  Once  the  corrections  were  made,  Pfeiffer’s  prediction  of  the  droplet 
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Figure  3.2  Comparison  of  Clewell’s  and  Pfeiffer’s  calculations  of  JP-4  droplet  tem¬ 
perature  versus  time  showing  the  error  in  Pfeiffer’s  calculation  of  tem¬ 
perature.  Jettison  velocity  =  175m/s,  Too  =  0®  (7  at  the  ground, 
Dq  =  500  nm,  jettison  altitude  =  1500  m. 

temperature  is  smooth  and  more  closely  matching  that  of  Clewell.  Pfeiffer’s  corrected 
and  uncorrected  droplet  temperature  predictions  are  compared  to  Clewell’s  in  figure 
3.3.  Notice  that  after  the  corrections  are  made,  the  prediction’s  of  Pfeiffer  and 
Clewell  are  almost  identical.  The  remaining  difference  in  the  initial  temperature 
decrease  is  due  to  differences  in  the  time  step  adjustment  which  is  explained  below. 

The  predictions  of  Clewell  and  Pfeiffer  are  still  not  exactly  the  same.  This  may 
be  because  Pfeiffer  uses  a  different  approach  to  adjusting  the  time  step.  If  the  mass 
evaporated  during  any  time  step  is  greater  than  some  predefined  maximum,  which 
is  nominally  1%  of  the  initial  droplet  mass,  Clewell  simply  halves  the  computed 
changes  in  time,  altitude,  and  mass  until  the  amount  of  mass  evaporated  is  below 
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Figure  3.3  Comparison  of  Pfeiffer’s  corrected  and  uncorrected  calculations  of  JP-4 
droplet  temperature  versus  time  to  those  of  Clewell  showing  the  error 
in  Pfeiffer’s  original  calculation  of  temperature.  Jettison  velocity  = 
175m/s,  Too  =  0"  C  at  the  ground,  Dq  =  500  //m,  jettison  altitude 
=  1500  m. 


the  defined  maximum  amount.  Also,  the  temperature  change  is  not  halved  if  the 
steady-state  temperature  is  calculated  and  results  in  a  droplet  temperature  that  is 
closer  to  the  ambient  temperature  than  it  was  on  the  previous  time  step.  Pfeiffer 
calculates  the  factor  needed  to  reduce  the  change  in  mass  evaporated  to  the  maximum 
amount  allowed  and  then  scales  the  changes  in  mass,  altitude,  temperature,  and  time 
uniformly  rather  than  just  halving  the  changes.  Furthermore,  Pfeifier  scales  the 
change  in  temperature  regardless  of  the  result  of  the  temperature  change.  Pfeiffer’s 
method  is  more  efficient  because  Clewell’s  method  may  require  more  than  one  scaling 
iteration  and  may  reduce  the  changes  in  mass,  temperature,  and  altitude  by  more 
than  is  necessary.  Figure  3.4  is  a  comparison  of  Clewell’s  prediction  of  JP-4  droplet 


3-4 


temperature,  using  Pfeiffer’s  method  of  scaling  the  changes,  to  that  of  Pfeiffer  with 
the  errors  corrected.  The  time  scale  has  again  been  shortened  to  the  first  10  seconds 
of  the  simulation  in  order  to  observe  the  point  where  the  largest  discrepancy  between 
the  two  models’  predictions  occurred.  Notice  that  Clewell’s  prediction  is  much  closer 
to  Pfeiffer’s  now  that  Pfeiffer’s  scaling  method  is  used  by  both  models. 
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Figure  3.4  Comparison  of  Pfeiffer’s  corrected  calculations  and  Clewell’s  calcula¬ 
tions  of  droplet  temperature  versus  time  using  Pfeiffer’s  method  of 
scaling.  Jettison  velocity  =  175m/s,  Too  =  0°  C  at  the  ground, 
Dq  =  500  fj,m,  jettison  altitude  =  1500  m. 


The  corrected  droplet  diameter  predictions  and  uncorrected  predictions  of 
Pfeiffer  are  compared  to  Clewell  predictions  of  droplet  diameter  using  Pfeiffer’s 
method  of  scaling  in  figure  3.5.  The  time  scale  has  been  shortened  to  the  first  200 
seconds  of  the  simulation,  where  the  temperature  difference  was  the  largest.  Pfeiffer’s 
corrected  temperature  calculation  led  to  a  lower  temperature  at  small  time  which 
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leads  to  a  larger  droplet  diameter.  After  Pfeiffer’s  code  is  corrected  and  Clewell’s 
code  is  modified  to  adjust  the  changes  as  Pfeiffer  does,  the  predictions  are  closer. 


Figure  3.5  Comparison  of  Pfeiffer’s  uncorrected  and  corrected  calculations  of 
droplet  diameter  versus  time  to  Clewell’s  calculation  using  Pfeiffer’s 
method  of  scaling.  Jettison  velocity  =  175m/s,  Too  =  0°  C  at  the 
ground,  Dq  =  500  ixm,  jettison  altitude  =  1500  m. 


3.3  Comparison  of  Clewell  and  Range,  Teske,  Polymeropolous 

For  comparison  to  experimental  data  found  in  (16),  where  the  droplets  were 
suspended  in  a  chamber  from  thin  thermocouple  wires  in  a  uniform  stream  of  dry 
air,  Clewell’s  algorithm  was  modified  to  match  the  experimental  conditions.  The 
algorithm  originally  computed  the  change  in  droplet  altitude,  the  initial  droplet  tem¬ 
perature  was  based  on  the  jettison  velocity,  and  the  Reynolds  number  was  computed 
using  the  droplets  terminal  velocity  at  its  current  size.  For  comparison  to  the  exper- 
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imental  measurements,  the  altitude  is  fixed  at  0  m  to  simulate  droplet  suspension 
and  the  Reynolds  number  is  calculated  using  a  fixed  air  velocity.  The  initial  droplet 
temperature  was  set  equal  to  the  ambient  temperature  instead  of  being  calculated. 

Figure  3.6  compares  the  measurements  of  heptane  evaporation  taken  by  Runge 
with  the  predictions  shown  in  Runge’s  thesis  (17)  and  Clewell  predictions.  The  initial 
conditions  are  reported  in  Table  3.1.  For  this  single  component  fuel,  the  pV  product 
was  assumed  by  Runge  to  be  constant  and  was  computed  to  be  8.5  x  10“®  g/(cni 
s)  by  using  the  Fuller  correlation  as  recommended  by  Reid  et  al.  (15).  The  Fuller 
correlation  was  also  used  to  compute  the  diffusivity  for  Clewell’s  prediction  as  well. 
For  both  predictions,  the  Lee-Kessler  relation  and  Pitzer  correlation  were  used  for 
the  calculation  of  the  vapor  pressure  and  latent  heat  of  vaporization,  respectively 
(15).  The  boiling  point,  molecular  weight,  and  density  of  heptane  were  taken  from 
the  data  corresponding  to  C7-parafl5ns  in  the  table  presented  in  the  appendix  of 
reference  (12). 


Constants 

Values 

Do{pm) 

638 

Tooi^C) 

20 

u(m/s) 

2 

pV{g/cmls) 

8.5x10-^ 

Table  3.1  Initial  conditions  for  evaporation  of  heptane. 

Figure  3.6  shows  that  the  predictions  of  Runge  are  almost  an  exact  match 
to  the  measured  data,  and  that  Clewell’s  predictions  are  reasonably  close.  This 
is  not  surprising  since  the  Runge’s  algorithm  was  developed  by  fitting  the  free  pT> 
parameter  to  JP-4  and  JP-8  evaporation  data  and  heptane  behaves  similarly. 

The  remaining  predictions  and  measurements  of  evaporation  that  appear  in 
Runge’s  thesis  cannot  be  compared  because  the  initial  diameters  of  the  droplets  are 
not  reported.  Runge  states  that  this  does  not  need  to  be  considered  because  the  data 
is  normalized  to  take  into  account  the  initial  droplet  size.  However,  this  is  not  the 
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Figure  3.6  Predictions  of  Heptane  evaporation  compared  to  experimental  data. 

Runge’s  prediction  and  the  experimental  data  were  taken  from  Runge’s 
thesis  and  Clewell’s  prediction  was  generated  using  Clewell’s  algorithm. 
Do  =  638iJ,m,  Too  =  20®  C,  u  =  2m fs 

case  unless  a  common  reference  Reynolds  number  is  used  for  all  calculations  because 
the  droplet  diameter  is  also  needed  to  compute  the  Reynolds  number. 

Unfortunately,  the  calculations  presented  in  reference  (16)  could  not  be  repro¬ 
duced  from  the  equations  in  the  appendix  of  reference  (16).  Recall  the  spherically 
symmetric  mass  vaporization  rate,  A4,  without  corrections  for  convective  effects,  as 
given  by  reference  (16): 

iVc 

M  =  -^R  pD  ln{l  -k  B) 

The  question  is,  what  is  pD?  There  is  no  p  or  D  listed  in  the  nomenclature  table 
of  reference  (16).  If  we  follow  the  notation  used  for  R,  then  D  =  D  '/D  q  and 
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p=z  p'Jp'^  where  the  primes  denote  dimensional  variables.  This  interpretation  does 
not  reproduce  the  predictions  though.  Alternatively,  pD  could  be  a  reduced  mass 
diffusivity-density  product,  {p  "D  ')/{p  'V  ')r,  since  the  variables  in  question  are  not 
primed  which  suggests  that  they  have  been  divided  by  a  reference  value.  However, 
this  approach  leads  to  a  dependence  on  the  value  of  {p  'V  ')r,  and  this  value  is  not 
specified.  Many  different  permutations  of  what  the  variables  could  be  were  tried  and 
none  led  to  predictions  similar  to  the  ones  presented  in  the  figures  of  the  article. 
Due  to  these  difficulties  with  the  equations  presented  in  the  article  (16),  the  curves 
presented  in  the  figures  are  used  rather  than  reproducing  them  using  the  equations. 


(tv/Ro^)x10‘^ 


Figure  3.7  Comparison  of  the  predictions  of  Clewell,  Runge’s  thesis,  and  Runge 
et  al.  to  the  experimentally  measured  non-dimensional  surface  area  of 
a  JP-4  droplet  versus  non-dimensional  time.  Ambient  air  velocity  =  3 
m/s.  Too  =  21"  C,  Do  =  636  pm. 
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Figures  3.7  through  3.14  are  comparisons  of  the  calculations  performed  using 
the  equations  of  Clewell,  Range’s  thesis,  and  the  Runge  et  al.  article  to  experimental 
data  measured  by  Runge  et  al.  The  non-dimensional  time  is  computed  as  10“^  tvjrl 
where  v  is  the  kinematic  viscosity  of  air.  The  non-dimensional  surface  area  is  com¬ 
puted  as  {D/DqY.  The  ambient  air  flow  is  denoted  by  «,  the  ambient  temperature 
is  denoted  by  Too,  and  the  initial  droplet  diameter  is  denoted  by  Dq. 

Notice  that  the  equations  used  in  the  thesis  by  Runge  (17)  produce  results  that 
are  very  close  to  those  presented  in  the  article.  Although  the  formulations  are  not 
the  same,  both  sets  of  equations  are  based  on  the  approach  of  Law  and  we  would 
expect  that  the  calculations  would  be  similar. 


Figure  3.8  Comparison  of  the  predictions  of  Clewell,  Range’s  thesis,  and  Runge 
et  al.  to  the  experimentally  measured  non-dimensional  surface  area  of 
a  JP-4  droplet  versus  non-dimensional  time.  Ambient  air  velocity  =  3 
m/s.  Too  =  —15'’  C,  Dq  =  600  fim. 
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Figure  3.9  Comparison  of  the  predictions  of  Clewell,  Runge’s  thesis,  and  Runge 
et  al.  to  the  experimentally  measured  non-dimensional  surface  area  of 
a  JP-8  droplet  versus  non-dimensional  time.  Ambient  air  velocity  =  3 
m/s,  Too  =  —14.5°  C,  Dq  =  645  urn. 

With  the  exception  of  figure  3.9,  Clewell  predicts  a  smaller  droplet  surface 
area  than  was  measured  and  calculated  by  the  other  approach.  The  main  difference 
between  the  two  approaches  is  the  calculation  of  the  gas  phase  diffusivity  and  mass 
density  of  the  fuel  components. 

Clewell  computes  the  diffusivity  of  each  species  in  air  using  an  empirical  cor¬ 
relation  that  is  derived  from  one  found  in  Bird  et  al.  (1)  using  parameters  that  are 
typical  of  hydrocarbons.  Clewell  handles  the  calculation  of  the  mass  density  as  in 
section  2.2.  Recall  that  the  gas  phase  mass  density  is  combined  with  the  mass  frac¬ 
tion  of  species  i  to  yield  the  gas  phase  mass  density  of  species  i.  The  mass  density 
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of  species  i  can  then  be  computed  by: 


pi  =  cxm 

The  approach  of  Runge  and  Runge  et  al.  requires  that  the  coefficients  used  to 
calculate  the  mass  density  times  diffusivity  product,  {pV),  be  determined  through 
fitting  curves  to  the  experimental  data.  A  reference  value,  {p'D)r,  is  used  for  the 
product  in  the  calculations  of  the  rates  of  change.  Then  the  dimensional  time  is 
scaled  by  the  factor  {p'D)r/{pT>)  in  order  to  match  the  experimental  data. 

This  difference  explains  why  the  calculations  of  Clewell’s  algorithm  appear 
to  be  more  sensitive  to  temperature  changes.  Clewell  computes  the  diffusivity  as  a 
function  of  temperature,  but,  Runge’s  algorithm  uses  the  same  pV  product  regardless 
of  the  temperature.  This  explains  why  Clewell  predicts  a  larger  JP-8  droplet  surface 
area  in  figure  3.9,  where  the  ambient  temperature  is  —14.5®  C,  and  a  smaller  JP-8 
droplet  surface  area  at  the  higher  ambient  temperatures.  The  same  effect  for  JP-4 
evaporation  is  seen  by  comparing  figure  3.7,  where  Too  =  21®  C,  and  figure  3.8, 
where  Too  =  —15®  C.  The  difference  between  the  prediction  of  Clewell  and  the 
measured  data  in  figure  3.8  is  greater  than  in  figure  3.7  because  the  diffusivity  is 
much  smaller  at  the  lower  temperature,  which  leads  to  a  larger  droplet  surface  area. 

In  each  comparison  of  the  measurements  to  the  calculations,  the  approach 
of  Runge  et  al.  and  Runge’s  thesis  produces  very  accurate  predictions  of  the  non- 
dimensional  surface  area.  This  is  to  be  expected  since  Runge  et  al.  used  experimental 
data  to  fit  the  parameters  used  in  the  calculation. 

Clewell’s  approach  was  not  adjusted  to  match  the  experimental  data  and  does 
not  produce  results  that  are  as  accurate  as  the  other  algorithms  that  use  a  pT> 
product  that  has  been  fit  to  experimental  data.  Although  the  number  of  comparisons 
to  experimental  data  is  small,  Clewell’s  algorithm  seems  to  produce  more  accurate 
predictions  of  JP-8  evaporation  than  JP-4  evaporation.  This  may  be  because  JP-8 
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Figure  3.10  Comparison  of  the  predictions  of  Clewell,  Runge’s  thesis,  and  Runge 
et  al.  to  the  experimentally  measured  non-dimensional  surface  area  of 
a  JP-8  droplet  versus  non-dimensional  time.  Ambient  air  velocity  = 
3  m/s,  Too  =  21°  C,  Do  =  636  fj,m. 

is  less  volatile  than  JP-4  and  the  JP-8  droplet  temperature  is  not  changed  through 
evaporative  cooling  as  much  as  the  JP-4  droplet  temperature  (compare  figures  3.13 
and  3.14). 

Figure  3.13  is  a  comparison  of  the  predictions  and  experimental  measurements 
of  JP-8  droplet  temperature  during  evaporation.  The  predictions  of  Clewell  and 
Runge  et.  al  appear  to  match  fairly  well  and  they  are  producing  the  same  qualitative 
behavior  as  the  measurements  indicate.  However,  the  experimental  data  shows  a 
slower  decrease  in  temperature  than  predicted  by  the  models.  This  may  be  caused 
by  the  effects  of  suspending  the  droplets  from  a  wire  because  the  wire  is  being  heated 
and  cooled  with  the  droplet.  This  increase  the  mass  being  cooled  by  the  evaporative 
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Figure  3.11  Comparison  of  the  predictions  of  Clewell,  Runge’s  thesis,  and  Runge 
et  al.  to  the  experimentally  measured  non-dimensional  surface  area  of 
a  JP-8  droplet  versus  non-dimensional  time.  Ambient  air  velocity  = 
1  m/s,  Tcx)  =  10°  C,  Do  =  647  iim. 

cooling  and  slow  the  process.  The  article  of  Runge  et  al.  (16)  examined  the  effect  of 
the  support  wire  on  droplet  evaporation  measurements.  The  support  wire  alters  the 
rate  of  heat  transfer  to  the  droplet  and  affects  the  droplet  internal  liquid  motion. 
They  performed  a  series  of  droplet  evaporation  measurements  using  different  support 
materials  and  sizes  until  there  was  no  apparent  effect  on  the  measured  results.  They 
concluded  that  a  25  fim  diameter  type  K  thermocouple  support  wire  had  negligible 
effects  on  the  droplet  evaporation.  However,  they  did  not  specify  whether  they 
considered  the  effect  on  the  droplet  temperature.  It  appears  as  though  the  wire  does 
affect  the  droplet  temperature  change  while  not  significantly  affecting  the  droplet 
diameter. 
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Figure  3.12  Comparison  of  the  predictions  of  Clewell,  Runge’s  thesis,  and  Runge 
et  al.  to  the  experimentally  measured  non-dimensional  surface  area  of 
a  JP-8  droplet  versus  non-dimensional  time.  Ambient  air  velocity  = 
3  m/s,  Too  =  10°  C,  Do  =  639  fim. 

Notice  that  the  scale  of  figure  3.13  is  only  1.2  K.  This  is  because  JP-8  is  less 
volatile  than  JP-4  and  so  the  temperature  is  not  affected  by  evaporative  cooling  as 
strongly.  This  means  that  although  the  temperature  predictions  for  JP-8  are  not 
exactly  what  was  measured,  the  predictions  are  only  wrong  by  a  small  amount. 

Figure  3.14  is  a  comparison  of  the  predictions  and  experimental  measurements 
of  JP-4  droplet  temperature  during  evaporation.  The  initial  temperature  decrease 
predicted  by  the  models  occurs  more  quickly  and  is  larger  than  the  experimental 
measurement.  It  appears  as  though  the  support  wire  has  again  slowed  and  damped 
the  initial  decrease  in  temperature.  Clewell’s  temperature  calculations  are  closer  to 
the  measured  values  than  the  calculations  of  Runge  et  al. 
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Figure  3.13  JP-8  droplet  temperature  predictions  of  Clewell  and  Runge  et  al.  com¬ 
pared  to  measured  JP-8  temperature  during  evaporation.  Too  =  O'*  C, 
Do  =  639  ixm,  u  =  3m/ s. 


Figures  3.15  through  3.18  are  comparisons  of  the  calculations  of  Clewell  and 
Runge’s  thesis  equations  to  experimental  measurements  of  JP-4  droplet  evaporation 
rates  measured  by  the  Arnold  Engineering  Development  Center  and  presented  in 
reference  (4).  In  each  case,  the  relationship  between  the  calculations  of  the  two 
approaches  is  the  same.  Clewell’s  algorithm  calculates  evaporation  rates  that  are 
in  agreement  with  than  those  calculated  by  Runge’s  thesis  algorithm  at  a  point 
in  time  between  20  and  30  seconds  into  the  simulations.  Before  the  time  when 
the  algorithms  agree,  Clewell’s  algorithm  calculates  higher  evaporation  rates  than 
Runge’s  algorithm.  At  times  beyond  the  point  of  agreement,  Clewell’s  algorithm 
calculates  lower  evaporation  rates  than  Runge’s  algorithm. 
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Figure  3.14  JP-4  droplet  temperature  predictions  of  Clewell  and  Runge  et  al.  com¬ 
pared  to  measured  JP-4  temperature  during  evaporation.  Too  =  0"  C, 
Do  =  646//m,M  =  Sm/s. 
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Figure  3.15  Comparison  of  calculations  of  Clewell  and  Runge’s  thesis  equations 
compared  to  measured  JP-4  evaporation  rate.  Too  =  20°  C,  Dq  = 
1235/im,  u  =  3m/ s. 
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Figure  3.16  Comparison  of  calculations  of  Clewell  and  Runge’s  thesis  equations 
compared  to  measured  JP-4  evaporation  rate.  Too  =  20°  C,  Do  = 
1347^7n,  u  =  Zm/s. 
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Figure  3.17  Comparison  of  calculations  of  Clewell  and  Runge’s  thesis  equations 
compared  to  measured  JP-4  evaporation  rate.  Too  =  20"  C,  Dq  = 
1060fj,m,u  =  3m/ s. 
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Figure  3.18  Comparison  of  calculations  of  Clewell  and  Runge’s  thesis  equations 
compared  to  measured  JP-4  evaporation  rate.  Too  =  20°  C,  Dq  — 
1179/im,u  =  Sm/s. 
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IV.  Surface  Evaporation 


4.1  Overview 

We  consider  the  evaporation  of  fuel  from  a  flat  surface  because,  in  toxicology 
studies,  this  flat  surface  could  be  the  skin  of  ground  crew  exposed  to  the  fuel. 

4.2  Assumptions 

•  The  surface  from  which  the  jet  fuel  is  evaporating  is  flat. 

•  The  surface  is  large  enough  that  we  only  have  to  consider  transport  normal  to 
the  surface. 

•  The  fuel  does  not  penetrate  the  surface. 

•  The  evaporating  fuel  is  uniformly  distributed  in  composition  and  temperature. 

•  Each  species  evaporates  independently  of  the  rest  of  the  mixture. 

•  Air  is  insoluble  in  the  fuel. 

•  The  gas  phase  concentration  of  i  at  the  surface  of  the  fuel  is  the  equilibrium 
vapor  pressure  of  species  i  in  the  mixture. 

•  The  vapor  pressure  of  each  species  in  the  mixture  follows  Raoult’s  Law. 

•  The  mixture  of  air  and  fuel  gases  at  the  surface  of  the  liquid  fuel  form  an  ideal 
mixture. 

•  The  concentration  of  fuel  in  the  ambient  airflow  is  zero. 

Now  we  have  simplified  the  problem  in  order  to  arrive  at  a  close  approximation 
to  the  solution.  The  validity  of  some  of  these  assumptions  needs  to  be  considered. 

The  assumption  that  the  skin  can  be  represented  by  a  flat  surface  may  be 
improved  by  considering  another  shape,  such  as  a  cylinder.  This  consideration  should 
be  addressed  in  future  work. 
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As  a  first  approximation,  the  assumption  that  fuel  species  evaporate  indepen¬ 
dently  of  the  rest  of  the  mixture  is  made  as  in  Lowell  (10).  In  section  4.3,  we  shall 
remove  this  assumption  and  consider  the  effects  that  the  species  mass  transfer  rates 
have  on  each  other  species.  At  that  time,  another  assumption  is  made  that  the  effec¬ 
tive  diffusivity  is  equal  to  the  diffusivity  in  a  binary  mixture  with  air.  This  second 
assumption  is  justified  when  the  diffusivity  is  only  slightly  position  dependent  (see 
reference  (1)). 

As  described  in  the  introduction,  fuel  is  coming  into  contact  with  the  ground 
crew.  Assuming  the  concentration  of  fuel  in  the  ambient  airflow  is  zero  may  not 
be  appropriate  during  the  initial  contact  period.  Therefore,  this  assumption  is  only 
valid  after  the  addition  of  fuel  to  the  skin  has  ceased. 

4.3  Equations 

Consider  the  heat  transfer  between  the  evaporating  liquid  fuel  and  the  air.  The 
heat  is  primarily  transferred  between  the  fuel  and  the  air  by  the  mechanisms  of  heat 
convection,  Qc,  and  radiation,  Qr-  The  convective  heat  flux  is  given  as  in  equation 
(2.27)  of  section  2.2: 

Qc  =  hAT  =  h{Too  —  Tp) 

where  h  is  the  heat  transfer  coefficient,  Tp  is  the  temperature  of  the  fuel,  and  T^o  is 
the  temperature  of  the  approaching  airflow.  The  heat  transfer  coefficient  is  given  as 
equation  (2.26)  in  section  2.2: 

where  k  is  the  thermal  conductivity  and  Ny,  is  the  Nusselt  number.  The  characteristic 
length  scale,  L,  is  the  length  of  the  surface  along  the  direction  of  the  ambient  airflow. 

The  Nusselt  number  is  a  function  of  the  Reynolds  and  Prandtl  numbers.  The 
coefficients  are  presented  in  reference  (2)  and  were  obtained  through  parameter  fit- 
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ting  experimental  data: 

Nu  =  OMN^^Nl/^ 

The  radiative  heat  flux  is  given  by  equation  (2.31)  in  section  2.2: 


Qr  =  cr^apToo  —  cfTf) 


where  a  is  the  Stefan-Boltzmann  constant,  ap  is  the  absorptivity  of  the  fuel,  and  ep 
is  the  emmisivity  of  the  fuel. 

Heat  is  also  lost  through  evaporative  cooling,  denoted  qh,  as  given  in  equation 
(2.32)  in  section  2.2: 

dm 

where  dm/dt  is  the  change  in  mass  with  time,  which  will  be  a  negative  number,  and 
AH  pap  is  the  latent  heat  of  vaporization.  The  product  is  divided  by  the  area  to  give 
the  evaporative  heat  loss  per  unit  area  per  unit  time.  The  total  heat  flux  is  the  sum 
of  these  contributions: 

q  =  Qc  +  Qr  +  qh  (4.1) 

The  change  in  temperature  with  respect  to  time  is  given  as  in  equation  (2.40): 

dTp  Aq 
dt  mCp 

where  A  is  the  surface  area  of  the  evaporating  fuel,  m  is  the  total  mass  and  Cp  is 
the  specific  heat  of  the  fuel. 

Consider  the  mass  transfer  of  a  pure  liquid.  A,  evaporating  from  a  flat  surface 
into  a  moving  gas  above  the  surface,  B.  The  mass  flux  of  gas  A  at  the  surface  of 
liquid  A  is  given  by  equation  (2.10)  in  section  2.2: 


nAo  =  YAoiriAo  +  nao)  +  ^AYa 


(4.3) 
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where  the  definitions  are  the  same  as  in  section  2.2.  Next,  we  follow  the  develop¬ 
ment  in  section  2.2  with  one  exception.  The  mass  flux  with  respect  to  stationary 
coordinates  is  used  so  that  the  diffusion  induced  bulk  flow  term  remains.  Then  sub¬ 
scripting  this  equations  as  in  section  2.3.2,  the  time  rate  of  change  in  mass  of  species 
i  is  given  by: 


_  AkiYiQ 

^  ■ 

where  A  is  the  surface  area.  Notice  that  and  ki  are  dependent  on  the  fuel 
temperature. 

The  mass  transfer  coefficient  is  given  by  equation  (2.13),  k  =  pViaNs/L  where 
L  is  the  length  of  the  surface  along  the  direction  of  the  ambient  airflow,  and  Via 
is  the  diffusivity  of  species  i  in  air.  The  Sherwood  number  is  an  empirical  fit  to 
experimental  measurement  of  mass  transfer  from  a  flat  plate  in  a  flowing  liquid  (see 
reference  (2)).  As  was  stated  in  section  2.2,  it  is  a  function  of  the  Schmidt  and 
Reynolds  numbers: 

Ns  = 


Note  that  this  equation  for  laminar  flow  over  a  flat  surface  differs  from  the  equation 
(2.19),  which  is  for  laminar  flow  around  a  sphere. 

The  mass  fraction  of  species  i  at  the  surface  can  be  calculated  as  follows: 

_  XioWi 

^iO  —  — 

Ej=l  XjoWj  +  XaOWa 

where  Xoo  and  Wa  are  the  mole  fraction  at  the  surface  and  molecular  weight  of  air, 
respectively.  Then,  as  we  have  assumed,  the  mole  fraction  of  i  at  the  surface  is  equal 
to  the  equilibrium  concentration: 


Xio  = 


_  tyvap 


Po 


4-4 


where  Pi,vap  is  the  vapor  pressure  of  species  i  in  the  mixture  and  Pq  is  the  total 
pressure  at  the  surface.  The  vapor  pressure  of  species  i  is  calculated  using  Raoult’s 
Law,  equation  (2.35): 

P  = 

tyvap  A.1^  tyvap 

where  the  superscripts  zero  denotes  pure  i  and  the  i  denotes  the  liquid  phase. 

The  assumption  that  each  species’  evaporation  rate  is  independent  of  the  rest  of 
the  species’  evaporation  rates  can  be  discarded.  Then,  as  presented  in  reference  (1), 
the  multicomponent  mass  flux  of  species  i  due  to  ordinary  diffusion  is  conveniently 
given  in  terms  of  the  multicomponent  diffiisivity  of  i  in  an  ideal  mixture  with  the 
other  N  species: 

N 

Hi  =  -pVimVYi  +  Yi^  Hj  .  (4.5) 

i=i 

Since  we  are  still  assuming  that  the  surface  is  long  enough  that  diffusion  in  the  y-z 
plane,  parallel  to  the  surface,  is  negligible,  the  vector  equation  (4.5)  simplifies  to  a 
one-dimensional  equation: 

N 

fit  =  p'Diffi^Yi  -|-  Yi  Tij 
3=1 

where  Ay  is  the  composition  difference  in  the  direction  normal  to  the  surface.  Eval¬ 
uating  this  equation  at  the  surface  and  denoting  this  with  the  subscript  zero,  the 
equation  for  the  mass  flux  at  the  surface  of  the  evaporating  fuel  is: 


N 

fliO  ~  pPim^^i  d"  Ytfi  fljo 

3=1 


Then,  multiplying  through  by  the  surface  area  of  the  fuel  yields  the  evaporation  rate 
of  species  i: 


drrii  _  ^  druj 

dt  dt 

3=1 


■A.p'Dim^^i 


(4.6) 


where  the  area  is  the  length  times  the  width  of  the  fuel  on  the  surface. 
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Rearranging  this  equation  to  get  the  derivative  terms  on  the  left-hand  side 
results  in: 

Bringing  drrii/dt  out  of  the  summation  and  substituting  Alf  =  li,o  —  ii,oo  =  ^,o 
leads  to: 

S  ^  =  -ApVi^Yi., 

Dividing  through  by  lio  yields: 

1  —  drrii  ^  drrij  _  .  ^ 

~v  IT  ^  ~dr  ~ 

^i,o  dt  dt 

This  can  be  set  up  and  solved  as  a  matrix  equation  of  the  form: 


B^  =  pV 
at 


(4.7) 


where  B  is  the  matrix  of  the  coefficients  is: 


B  = 


f  (1  -  yi,o)/yi.o  1 

1  (1  -  y2.o)/l2,o 


V 


1 


(1  -  yiv,o)/yv,o ) 


Equations  (4.4)  and  (4.7)  have  been  implemented  for  comparison  to  the  exper¬ 
imental  data  below. 

To  describe  the  mass  flux  due  to  evaporation  of  fuel  from  a  surface,  equations 
(4.4)  or  (4.7)  are  appropriate.  The  mass  flux  with  respect  to  mass  average  velocity, 
ji,  should  not  be  used  to  describe  the  evaporation  from  a  flat  surface.  The  mass 
average  velocity  in  the  surface  evaporation  system  would  be  moving  away  from  the 
surface  since  it  has  been  assumed  that  no  mass  is  moving  through  the  surface.  As  we . 
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are  interested  in  the  mass  flux  at  the  surface  and  do  not  want  to  obtain  a  solution 
away  from  the  surface,  we  must  use  rij. 

When  applied  to  the  problem  of  evaporation  of  fuel  from  skin,  the  assumption 
that  no  mass  is  transported  through  the  surface  is  not  valid.  Dermal  absorption  does 
occur  and  this  problem  may  be  handled  by  also  computing  the  dermal  absorption 
rate  of  the  fuel  through  the  skin. 

Comparison  to  Experimental  Data 

Evaporation  of  Arctic  Diesel  40,  a  No.  2  fuel  oil,  has  been  measured  by  Regnier 
and  Scott  (14).  A  dark  evaporation  chamber  provided  a  constant  ambient  airflow  of 
21  km! hr  and  allowed  the  ambient  temperature  to  be  adjusted  to  5,  10,  20,  and  30 
degrees  C.  Note  that  the  high  ambient  flow  velocity  used  in  this  experiment  is  more 
characteristic  of  the  ambient  flow  velocity  encountered  during  a  cold  start  situation. 
The  vaporized  fuel  was  removed  from  the  air  stream  by  using  a  water  aspirator  at 
the  exhaust  outlet.  Samples  of  oil  weighing  12  ±  7  grams  were  measured  into  90 
mm  diameter  Petri  dishes  and  then  placed  in  the  chamber  and  weighed  periodically. 
The  oil  depth  was  approximately  3  mm  for  each  case. 

We  wish  to  determine  if  the  flat  surface  evaporation  model  predicts  the  evap¬ 
oration  of  fuel  from  a  flat  surface  more  accurately  than  the  droplet  evaporation 
models.  Therefore,  the  experimental  measurements  of  the  percentage  of  initial  mass 
of  oil  remaining  are  compared  to  the  flat  surface  evaporation  model  and  Clewell’s 
droplet  evaporation  model. 

For  the  flat  surface  evaporation  simulations,  the  surface  area  of  the  fuel  in  the 
circular  Petri  dish  was  computed  as  A  =  7rr^  =  7r(0.9)^,  cm^.  The  fuel  mass  was  set 
equal  to  12.0  g  to  match  the  experimental  conditions. 

For  the  spherical  droplet  evaporation  simulations,  the  diameter  was  calculated 
from  the  initial  droplet  mass,  which  was  set  equal  to  12.0  g.  The  volume  of  the 


4-7 


droplet  was  then  calculated  using: 


P 

where  m  is  the  mass  of  the  droplet  and  it  was  assumed  that  the  density  of  the  droplet, 
p,  was  equal  to  the  sum  of  the  densities  of  each  component  of  the  fuel.  The  volume 
of  a  spherical  droplet  is  given  by: 


which  was  solved  for  r  to  obtain  the  droplet  radius  required  to  result  in  the  given 
droplet  mass. 

The  fuel  oil  composition  from  reference  (14)  is  used  for  the  simulations  and  is 
shown  in  table  4.1.  The  species  listed  were  obtained  by  chromatography  and  only 
account  for  18%  of  the  mixture.  However,  it  is  hypothesized  by  Regnier  and  Scott 
that  the  composition  listed  in  table  4.1  is  sufficient  for  characterizing  the  evaporative 
behavior  of  the  oil.  This  follows  from  the  assumption  that  the  components  that  are 
lighter  than  nonane  will  evaporate  very  quickly  and  components  that  are  heavier 
than  octadecane  will  evaporate  very  little. 


Volume  Fraction 

Molecular  Weight 

Boiling  Point,  K 

Density  at  20°  C,  kg/rn^ 

0.039 

128.3 

415.15 

720.0 

BEH 

0.103 

142.3 

433.15 

720.0 

C-11 

156.3 

469.15 

740.0 

C-12 

170.3 

489.15 

750.0 

C-13 

0.155 

184.4 

760.0 

C-14 

0.158 

198.4 

527.15 

760.0 

C-15 

0.107 

212.4 

544.15 

770.0 

C-16 

0.073 

226.4 

560.15 

C-17 

0.044 

240.4 

575.20 

778.0 

C-18 

0.019 

254.4 

589.25 

777.0 

Table  4.1  Arctic  Diesel  40  characterization 
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Figures  4.1,  4.2,  4.3,  and  4.4  are  the  comparisons  of  the  predictions  of  the  flat 
surface  model  and  Clewell’s  droplet  evaporation  model  to  the  experimental  mea¬ 
surements  of  the  percent  of  initial  mass  remaining  at  30"  C,  20"  C,  10"  C,  and  5"  C, 
respectively.  The  thermodynamic  properties  used  for  the  predictions  are  obtained 
using  the  empirical  relations  for  hydrocarbons  that  are  derived  by  Clewell  in  reference 

(3).  . 

The  predictions  of  the  flat  surface  evaporation  model  are  in  reasonably  good 
agreement  with  the  measured  data.  For  the  higher  temperatures  though,  the  flat 
surface  model  predictions  of  the  mass  evaporated  are  too  high  at  later  times.  Al¬ 
though  the  initial  surface  area  and  mass  are  equivalent  for  the  droplet  and  the  fuel  in 
the  circular  Petri  dish,  the  droplet  evaporation  predictions  are  poor.  The  Sherwood 
number,  Ng,  is  computed  in  the  droplet  model  by: 


Ns  =  2.0  +  0.6N]/^Ny^ 

while  Ns  is  computed  by  the  flat  surface  model  by: 

Ns  =  OMN^/^N^^ 


Given  the  differences  in  the  calculations  of  the  Sherwood  numbers,  and  therefore  the 
mass  transfer  coefficients,  we  would  expect  the  droplet  model  to  predict  evaporation 
rates  that  are  too  large. 

At  large  times,  the  predictions  of  the  flat  surface  model  at  the  two  lower  tem¬ 
peratures  (figures  4.3  and  4.4)  are  in  closer  agreement  with  the  measurements  than 
the  predictions  at  the  higher  temperatures  (figures  4.1  and  4.2).  At  small  times,  the 
predictions  at  the  two  higher  ambient  temperatures  are  closer  to  the  experimental 
measurements  than  the  predictions  at  the  lower  ambient  temperatures.  When  the 
temperature  is  5  degrees  Celsius,  the  surface  evaporation  model  under-predicts  the 
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Figure  4.1  Comparison  of  predictions  of  percentage  of  initial  mass  of  oil  remaining 
to  experimentally  measured  values  at  T  =  30°  C. 

amount  of  fuel  evaporated  by  as  much  as  5%  in  the  mid-stages  of  the  simulation, 
and  then  comes  close  to  the  measured  mass  again  at  later  times. 

The  error  in  the  flat  surface  model’s  predictions  of  the  mass  remaining  may  be 
because  we  have  assumed  that  we  have  a  flat  plate  wetted  with  fuel.  However,  the 
experimental  data  is  the  measurement  of  the  evaporation  of  the  fuel  from  inside  of  a 
Petri  dish.  As  the  fuel  evaporates,  the  fuel  surface  falls  below  the  edges  of  the  dish, 
which  would  affect  the  flow  of  air  across  the  surface.  Secondly,  the  air  velocity  is  very 
high,  producing  a  Reynolds  number  on  the  order  of  3.0  x  10®  which  indicates  that 
we  are  in  the  transition  regime.  Turbulent  effects,  which  have  not  been  modeled, 
are  likely  affecting  the  evaporation  rate.  Another  possible  explanation  is  that  the 
predicted  fuel  temperature  is  higher  than  the  actual  temperature.  We  have  followed 
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Figure  4.2  Comparison  of  predictions  of  percentage  of  initial  mass  of  oil  remaining 
to  experimentally  measured  values  at  T  =  20®  C. 

the  approach  of  Clewell  which  uses  empirical  relations  for  hydrocarbons  that  air 
fairly  general.  However,  the  linear  formula  for  the  latent  heat  of  vaporization  was 
derived  by  Clewell  from  inspection  of  the  heats  of  vaporization  for  the  range  of 
hydrocarbons  that  are  in  JP-4.  JP-4  is  composed  of  several  hydrocarbons  that  are 
lighter  than  those  listed  in  table  4.1,  which  may  lead  to  a  relationship  for  the  latent 
heat  that  is  too  small.  This  could  lead  underestimating  the  amount  of  evaporative 
cooling,  producing  temperature  predictions  that  are  too  high.  In  addition,  the  higher 
temperatures  would  increase  the  error  in  Clewell’s  approach  for  predicting  the  latent 
heat  of  vaporization.  Without  measured  temperature  data,  this  cannot  be  confirmed 
though. 
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Time  (minutes) 

Figure  4.3  Comparison  of  predictions  and  measurements  of  percentage  of  initial 
mass  of  oil  remaining  at  T  =  10°  C. 

It  should  be  noted  here  that  Regnier  and  Scott  developed  an  activation  energy 
evaporation  model  (14).  The  predictions  of  their  model  were  more  accurate  than 
both  the  flat  surface  or  spherical  model  discussed  earlier.  This  suggests  that  an 
activation  energy  approach  to  calculating  the  evaporation  of  fuel  from  skin  should 
be  considered  in  future  work. 
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V.  Summary  and  Conclusions 


5.1  Summary 

For  this  research,  two  previously  developed  fuel  droplet  evaporation  models 
were  compared  and  the  equations  to  describe  surface  fuel  evaporation  were  devel¬ 
oped  in  order  to  predict  evaporation  of  fuel  from  skin  following  exposure.  The 
equations  of  the  two  previous  algorithms  were  shown  to  be  derived  from  basic  prin¬ 
ciples  of  heat  and  mass  transfer  and  the  predictions  were  compared  to  experimental 
laboratory  data.  The  surface  evaporation  equations  were  then  derived  from  those 
basic  principles,  following  the  approach  of  Clewell,  and  its  predictions  of  the  surface 
evaporation  was  compared  to  experimental  data.  This  research  resulted  in  a  com¬ 
puter  program  that  may  be  used  to  predict  the  evaporation  of  fuel  from  a  surface. 
The  surface  evaporation  algorithm  can  be  used  in  conjunction  with  a  program  that 
predicts  the  dermal  absorption  of  fuel  to  determine  the  dose  of  fuel  delivered  to  the 
exposed  ground  crew  workers. 

5.2  Conclusions 

The  two  previous  descriptions  produce  predictions  that  agree  very  well  with 
their  respective  sets  of  experimental  data.  However,  the  predictions  presented  in 
reference  (16)  could  not  be  reproduced  from  the  equations  shown  in  the  appendix 
of  reference  (16).  If  the  problem  with  the  inconsistent  notation  used  to  present  the 
equations  is  resolved  and  the  predictions  that  appear  in  the  article  can  be  reproduced, 
then  the  approach  of  Runge  et  al.  produce  more  accurate  predictions  of  JP-4  and 
JP-8  evaporation  than  the  approach  of  Clewell.  The  equations  used  in  Runge’s  thesis 
(reference  (17))  may  be  used  to  produce  results  that  match  those  of  reference  (16) 
almost  exactly. 

In  terms  of  adapting  to  meet  the  changing  needs  of  the  Air  Force  through  the 
modeling  of  new  aircraft  fuels,  the  approach  of  Clewell  is  superior.  The  only  exper- 
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iments  which  must  be  performed  to  model  a  new  hydrocarbon  fuel  using  Clewell’s 
approach  is  one  to  ascertain  the  predominant  species  with  which  to  characterize  the 
fuel’s  behavior.  It  has  already  been  partly  demonstrated  by  the  predictions  of  the 
evaporation  of  the  Arctic  Diesel  40  that  this  approach  is  robust  while  producing  rea¬ 
sonably  accurate  results.  Next,  the  fuel  specific  correlations  for  the  thermodynamic 
properties,  such  as  the  specific  heat  and  latent  heat,  could  be  recalculated  using  data 
that  can  be  found  in  the  literature.  To  model  the  evaporation  of  a  new  fuel  using  the 
approach  of  Runge  et  al.,  one  would  have  to  perform  several  evaporation  experiments 
in  order  to  obtain  constants  for  use  in  the  calculation  of  the  pV  product. 

The  calculations  of  the  algorithm  are  within  5%  of  the  measured  mass  remain¬ 
ing  for  temperatures  of  5  and  10  degrees  Celsius  using  a  10  species  characterization 
of  the  fuel  that  accounted  for  approximately  18%  of  the  actual  fuel  oil  mixture. 
For  higher  temperature  vaporization,  the  surface  evaporation  algorithm  reproduced 
the  experimental  data  to  within  10%  mass  remaining  for  temperatures  of  20  and  30 
degrees  Celsius. 

5.3  Recommendations  for  Further  Research 

Further  research  into  the  evaporation  of  jet  fuel  from  skin  is  needed.  It  remains 
to  be  determined  if  the  approximation  of  the  skin  as  a  flat  surface  is  a  good  assump¬ 
tion,  or  if  there  is  a  better  geometry  such  as  a  cylinder.  Secondly,  the  performance 
of  the  model  in  high  ambient  flow  velocities,  such  as  those  encountered  during  cold 
start,  has  been  investigated  in  chapter  IV.  However,  turbulent  effects  at  the  high 
ambient  flow  velocity  should  be  considered  in  future  work. 

The  activation  energy  approach  used  by  Regnier  and  Scott  (14)  produced  more 
accurate  predictions  of  the  mass  evaporated  than  the  other  models  compared.  Re¬ 
search  into  an  activation  energy  approach  to  calculating  the  evaporation  of  fuel  from 
skin  should  be  conducted. 
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